Runge-Kutta-Verfahren mit Schrittweitensteuerung

am Beispiel der DG

Initialisierung

> restart:

> with(plots):

Eingabedaten

( hier dürfen Veränderungen vorgenommen werden )

> x0:=0: Intervallanfang

> x1:=5: Intervallende

> alias(y=y(x), y0=y(x0)): Definition

> initval:= y0=-4: Anfangswert:=-4: Anfangswert

> h:=0.1: Schrittweite

> f:=(x,y) -> y^2*sin(x): zu lösende Funktion

Ausgangsdgl

> dgl:=diff(y,x)=f(x,y); Ausgabe der zu lösenden Differentialgleichung

Exakte Lösung

> exakt:=dsolve({dgl,initval},y); exakte Lösung der DGL

> y_e:=unapply(rhs(%),x):

Verfahren 4. Ordnung mit Schrittweitensteuerung

> xh[0]:= x0:
yh[0]:= evalf(Anfangswert):
ye[0]:= evalf(Anfangswert):
i:=0:

> while xh[i] < x1 do
i:=i+1:
do
K1:=evalf(f(xh[i-1],yh[i-1]));
K2:=evalf(f(xh[i-1]+h/2,evalf(yh[i-1]+K1*h/2))):
K3:=evalf(f(xh[i-1]+h/2,evalf(yh[i-1]+K2*h/2))):
kappa:=evalf(abs((K3-K2)/(K2-K1)));
if i = 1 then
kappalt := kappa:
end if:
if 0.5*(kappa + kappalt) < 0.05 then
h:= evalf(h * 2);
elif 0.5*(kappa + kappalt) > 0.2 then
h:= evalf(h / 2);
else
if (evalf(x1-xh[i-1]) < h) then h:= evalf(x1-xh[i-1]) end if:
break:
end if:
end do;
k[i]:=kappa:
kappalt:=kappa:
hs[i]:=h:
K4:=evalf(f(xh[i-1]+h,yh[i-1]+K3*h)):
xh[i]:=evalf(xh[i-1]+h):
yh[i]:=evalf(yh[i-1]+h/6*(K1+2*(K2+K3)+K4)):
ye[i]:=evalf(y_e(xh[i])):
AbsFehler[i]:=ye[i]-yh[i]:
end do:

> imax:=i:

> out[-1,1]:=k: out[-1,2]:=X[k]: out[-1,3]:=Y[k]:out[-1,4]:=Exakt: out[-1,5]:=AbsFehler: out[-1,6]:=Schrittkennzahl: out[-1,7]:=Schrittweite:

> for i from 0 by 1 to imax do Schleifen Anfang für Ausgabe der Lösungen

> out[i,1]:=i+1: out[i,2]:=xh[i+1]: out[i,3]:=yh[i+1]: out[i,4]:=ye[i+1]:out[i,5]:=AbsFehler[i+1]:out[i,6]:=k[i+1]:out[i,7]:=hs[i+1]:

> od: Schleifen Ende

Ausgabe

> werte:=false; Ausgabe der Werte (true=>Ausgabe ; false=>keine Ausgabe)

> if (werte=true) then a:=(i,j) -> out[i-2,j]: matrix((imax+1),7,a) else end if;

> x:='x': x wird zurückgesetzt

> plot_exakt:=plot(y_e(xs), xs=x0..x1,legend=[Exakt]): Ausgabe des exakten Plots

> plott:=[xh[1],yh[1]]:plotfehler:=[xh[1],AbsFehler[1]]:plotk:=[xh[1],k[1]]:

>

> for i from 2 to imax do

> plott:=plott,[xh[i],yh[i]]:
plotfehler:=plotfehler,[xh[i],AbsFehler[i]]:
plotk:=plotk,[xh[i],k[i]]

> od:

> plott:=[plott]:plotfehler:=[plotfehler]:plotk:=[plotk]:

> plot_num:=plot(plott,style=point,color=[blue],symbol=cross,legend=[Runge_Kutta]):
plot_err:=plot(plotfehler,style=point,color=[blue],symbol=circle,legend=[Fehler]):
plot_k:=plot(plotk,color=[blue],symbol=cross,legend=[Kappa]):

> display(plot_exakt,plot_num);
display(plot_err);
display(plot_k);