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

> restart;
with(plots):

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

f := proc (t, y) options operator, arrow; -2*y end proc

> tk := a + dt*k;

tk := a+dt*k

Exakte L?sung

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

y(t) = exp(-2*t)

> subs(t=1,%);

y(1) = exp(-2)

> ex := rhs(%);

ex := exp(-2)

> evalf(ex);

.1353352832

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

[Plot]

Euler

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

phi := proc (t, y, dt) options operator, arrow; f(t, y) end proc

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

y(k+1) = y(k)-2*dt*y(k)

> collect(%,y(k));

y(k+1) = (1-2*dt)*y(k)

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

(1-2*dt)^k

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

eu := proc (k, dt) options operator, arrow; (1-2*dt)^k end proc

> 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]);

[Plot]

>

Heun

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

phi := proc (t, y, dt) options operator, arrow; 1/2*f(t, y)+1/2*f(t+dt, y+dt*f(t, y)) end proc

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

y(k+1) = y(k)-2*dt*y(k)+2*dt^2*y(k)

> collect(%,y(k));

y(k+1) = (1-2*dt+2*dt^2)*y(k)

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

(1-2*dt+2*dt^2)^k

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

he := proc (k, dt) options operator, arrow; (1-2*dt+2*dt^2)^k end proc

> 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]);

[Plot]

>

Mittelpunktsregel

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

y(k+1) = y(k-1)-4*dt*y(k)

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

-1/2*(2*(1/(2*dt-(4*dt^2+1)^(1/2)))^k*dt-(1/(2*dt-(4*dt^2+1)^(1/2)))^k*(4*dt^2+1)^(1/2)+(1/(2*dt-(4*dt^2+1)^(1/2)))^k*exp(-2*dt)-2*(1/(2*dt+(4*dt^2+1)^(1/2)))^k*dt-(1/(2*dt+(4*dt^2+1)^(1/2)))^k*(4*dt^...

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

cp := x^2+4*dt*x-1

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

xi := [-2*dt+(4*dt^2+1)^(1/2), -2*dt-(4*dt^2+1)^(1/2)]

> xi[1];

-2*dt+(4*dt^2+1)^(1/2)

> xi[2];

-2*dt-(4*dt^2+1)^(1/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]);

[Plot]

>

Konvergenz Euler

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

l1 := [0., 0.6250000000e-1, .1001129150, .1180670870, .1267887864, .1310840325, .1332151629, .1342765997, .1348062856, .1350708705]

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

l2 := [-.1353352832, -0.7283528320e-1, -0.352223682e-1, -0.172681962e-1, -0.85464968e-2, -0.42512507e-2, -0.21201203e-2, -0.10586835e-2, -0.5289976e-3, -0.2644127e-3]

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

l3 := [1.858100597, 2.067870133, 2.039724809, 2.020499932, 2.010348813, 2.005193149, 2.002600683, 2.001301140, 2.000651255]

>

Konvergenz Heun

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

l1 := [.2500000000, .1525878906, .1387778781, .1361117386, .1355200935, .1353803906, .1353464272, .1353380529, .1353359736, .1353354556]

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

l2 := [.1146647168, 0.172526074e-1, 0.34425949e-2, 0.7764554e-3, 0.1848103e-3, 0.451074e-4, 0.111440e-4, 0.27697e-5, 0.6904e-6, 0.1724e-6]

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

l3 := [6.646225358, 5.011512508, 4.433731673, 4.201364318, 4.097117103, 4.047684853, 4.023540456, 4.011732329, 4.004640371]

>

>