################################################################################
# Worksheet: de-soln2
# CP 1. Pendulum behavior. 
# http://cadigweb.ew.usna.edu/~wdj/teach/sm212/sm212_extra_credit/sm212_extra_credit.html
################################################################################


%hide
%latex 
For a simple pendulum without an external driving force, the angle $\theta(t)$ 
\newline
the pendulum makes with the vertical is described by the DE 
\newline
$ \theta'' + {{g}\over{L}}\sin \theta = 0, $ and solve it using 
\newline 
Maxima's \verb+de_solve+ command.

%sage
de = lambda x: 'diff(theta,t,2) + (%s/%s)*sin(theta)=0%s+%s'%(x[0],x[1])
de([32.2,8])

%latex
Here x[0] represents the gravitational constant g and x[1] the length L.
Next, we input the 2nd order DE above as a system of 2 1st order DEs.

%sage
f = lambda z : z[2]
g = lambda z : -(32.2/8)*sin(z[1])

%This is now plotted, using the function \verb+eulers_method_2x2_plot+
\newline
entered below (because to avoid extra blank lines and wrong indentations
\newline
when pasting this is).

%sage
P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0)

%sage
P[0].show()

%hide
%latex 
Note that since the DE is non-linear, the usual de_solve command
\verb+maxima.de_solve(de([32.2,8]),["t",'theta'],[0,3/4,0])+
does note work.

########### need a numerical plot of the soln....... ###############

%hide
%latex
where $g = 32.2 {\rm ft/sec^2},$ and $L$ is the length of the pendulum in ft. 
\newline
This needs to be solved by numerical means to get a good approximation.

   1. Let $L = 8$, $\theta (0) = \pi/4$, $\theta '(0) = 0$. 
For $0 \leq t \leq 5\pi,$ 
\newline
plot a solution to this pendulum problem.

   2. If you have had any experience with a pendulum, does this solution match you experience?


################ This needs the following Python function. ################


def eulers_method_2x2_plot(f,g, t0, x0, y0, h, t1):
    """
    This plots the soln in the rectangle 
    (xrange[0],xrange[1]) x (yrange[0],yrange[1])
    and plots using Euler's method the numerical solution of 
    the 1st order ODEs x' = f(t,x,y), x(a)=x0, y' = g(t,x,y), y(a) = y0. 
    
    The following example plots the solution to theta''+sin(theta)=0,
    theta(0)=3/4, theta'(0) = 0.  Type P[0].show() to plot the solution, 
    (P[0]+P[1]).show() to plot (t,theta(t)) and (t,theta'(t)).
   
    EXAMPLES:
        sage: f = lambda z : z[2]; g = lambda z : -sin(z[1])
        sage: P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0)
    """
    #print "%10s %20s %25s %20s %20s"%("t", "x","h*f(t,x,y)","y", "h*g(t,x,y)")
    n=int((1.0)*(t1-t0)/h)
    t00 = t0; x00 = x0; y00 = y0
    #print RR(f([t00,x00,y00])), RR(g([t00,x00,y00]))
    soln = [[t00,x00,y00]]
    for i in range(n+1):
        #print "%10r %20r %25r %20r %20r"%(t00,x00,h*f([t00,x00,y00]),y00,h*g([t00,x00,y00]))
        x01 = x00 + h*RR(f([t00,x00,y00]))
        y00 = y00 + h*g([t00,x00,y00])
        x00 = x01
        t00 = t00 + h
        #print i,t00,x00,y00
        soln.append([t00,x00,y00])
    #print soln
    Q1 = line([[x[0],x[1]] for x in soln], rgbcolor=(1/4,1/8,3/4))
    Q2 = line([[x[0],x[2]] for x in soln], rgbcolor=(1/2,1/8,1/4))
    return [Q1,Q2]