Modellbildung und Simulation, SS 2009
Blatt 8: numerische L?sung von ODEs

 > restart; with(plots):

 > f := (t,y) -> -2*y;

 > tk := a + dt*k;

Exakte L?sung

 > dsolve({diff(y(t),t)=-2*y(t),y(0)=1},y(t));

 > subs(t=1,%);

 > ex := rhs(%);

 > evalf(ex);

 > pex := plot(exp(-2*t),t=0..1,thickness=3): display(pex);

Euler

 > phi := (t,y,dt) -> f(t,y);

 > y(k+1) = y(k) + dt*phi(tk,y(k),dt);

 > collect(%,y(k));

 > rsolve({%,y(0)=1},y(k));

 > eu := unapply(%,k,dt);

 > peu := plot([seq([1/16*k,eu(k,1/16)], k=0..16)],            style=point, symbolsize=20, symbol=diamond,            color=blue, thickness=3): display([pex, peu]);

 >

Heun

 > phi := (t,y,dt) -> (f(t,y)+f(t+dt,y+dt*f(t,y)))/2;

 > y(k+1) = simplify(y(k) + dt*phi(tk,y(k),dt));

 > collect(%,y(k));

 > rsolve({%,y(0)=1},y(k));

 > he := unapply(%,k,dt);

 > phe := plot([seq([1/4*k,he(k,1/4)], k=0..4)],            style=point, symbolsize=20, symbol=diamond,            color=blue, thickness=3): display([pex, phe]);

 >

Mittelpunktsregel

 > y(k+1) = y(k-1) + 2*dt*f(tk,y(k));

 > simplify(rsolve({%,y(0)=1, y(1)=exp(-2*dt)},y(k)));

 > cp := x^2 + 2*dt*2*x - 1;

 > xi := [solve(cp,x)];

 > xi[1];

 > xi[2];

 > mp := proc(k, dt)  option remember;  if k=0 then return 1;  elif k=1 then return exp(-2*dt);  else return mp(k-2,dt)-4*dt*mp(k-1,dt);  end if; end proc:

 > pmp := plot([seq([1/4*k,mp(k,1/4)], k=0..4)],            style=point, symbolsize=20, symbol=diamond,            color=blue, thickness=3): display([pex, pmp]);

 >

Konvergenz Euler

 > l1 := [seq(evalf(eu(2^n,2^(-n))),n=1..10)];

 > l2 := [seq(evalf(eu(2^n,2^(-n))-ex),n=1..10)];

 > l3 := [seq(l2[n]/l2[n+1], n=1..9)];

 >

Konvergenz Heun

 > l1 := [seq(evalf(he(2^n,2^(-n))),n=1..10)];

 > l2 := [seq(evalf(he(2^n,2^(-n))-ex),n=1..10)];

 > l3 := [seq(l2[n]/l2[n+1], n=1..9)];

 >

 >