> Euler and Improved Euler with Maple

Euler and Improved Euler with Maple 


> f1:=(x,v)->v;

                                f1 := (x,v) -> v
> f2:=(x,v)->-g/l*sin(x);

                                             g sin(x)
                            f2 := (x,v) -> - --------
                                                 l
> h:=.05;

                                    h := .05
> l:=10;

                                     l := 10
> g:=9.8;

                                    g := 9.8
> xold:=Pi/4;

                                 xold := 1/4 Pi
> vold:=0;

                                    vold := 0
> soln:=0,evalf(Pi/4);

                             soln := 0, .7853981635
> t:=0;

                                     t := 0
> 

                             soln := 0, .7853981635

                                     t := 0
> for j from 1 to 800 do
> xnew:=evalf(xold+h*f1(xold,vold)):
> vnew:=evalf(vold+h*f2(xold,vold)):
> t:=t+h:
> soln:=soln,t,xnew:
> xold:=xnew:
> vold:=vnew:
> od:
> soln:
> plot1:=plot([soln]):
> t:='t';

                                     t := t
> de1:=diff(x(t),t)=v(t);

                                      d
                             de1 := ---- x(t) = v(t)
                                     dt
> de2:=diff(v(t),t)=-g/l*sin(x(t));

                            d
                   de2 := ---- v(t) =  - .9800000000 sin(x(t))
                           dt
> ss:=dsolve({de1,de2,x(0)=Pi/4,v(0)=0},{x(t),v(t)},numeric);

ss := proc(rkf45_x) ... end

> ss(2);

          [t = 2, x(t) = -.2598078029448924, v(t) = -.7129459288917021]
> with(plots):
> 


> 


> plot2:=odeplot(ss,[t,x(t)],0..40,numpoints=100):
> display([plot1,plot2]);
> xold:=Pi/4;

                                 xold := 1/4 Pi
> vold:=0;

                                    vold := 0
> soln:=0,evalf(Pi/4);

                             soln := 0, .7853981635
> t:=0;

                                     t := 0
> 
> for j from 1 to 800 do
> e1:=xold+h*f1(xold,vold):
> e2:=vold+h*f2(xold,vold):
> xnew:=evalf(xold+h*(f1(xold,vold)+f1(e1,e2))/2):
> vnew:=evalf(vold+h*(f2(xold,vold)+f2(e1,e2))/2):
> t:=t+h:
> soln:=soln,t,xnew:
> xold:=xnew:
> vold:=vnew:
> od:
> plot3:=plot([soln]):
> display([plot1,plot3]);
>