CV7 - Metoda Runge-Kutta

Program RK4 pro výpočet řešení diferenciální rovnice y'=f(t,y) na intervalu ⟨(Subscript[t, 0]),(Subscript[t, 1])⟩ s počáteční podmínkou y((Subscript[t, 0]))=(Subscript[y, 0]) s krokem dt metodou Runge-Kutta.

RK4[f_,t0_,t1_, y0_, dt_,tisk_] :=
	Module[{ k1, k2, k3, k4 },
		 ta=t0;
		ya=y0;
		vys={{ta,ya}};
		If[tisk==1,{Print["t       y"];
		Print[ta,"      ",ya]}];	
		While[t1-ta>0.00000001,
		k1 = dt N[f[t,y]/.{t->ta,y -> ya} ];
		k2 = dt N[f[t,y]/.{t->ta+dt/2,y -> ya+k1/2}] ;
		k3 = dt N[f[t,y]/.{t->ta+dt/2,y -> ya+k2/2}] ;
		k4 = dt N[f[t,y]/.{t->ta+dt,y -> ya+k3}];
		ya=N[ya + (k1 + 2 k2 + 2 k3 + k4)/6];
		ta=ta+dt;
		If[tisk==1,Print[ta,"      ",ya]];	
		AppendTo[vys,{ta,N[ya,16]}];	
	];
	vys
	]

Použítí programu RK4 a NDSolve pro výpočet řešení diferenciální rovnice

Příklad 1:
Mějme rovnici y'=t y na intervalu ⟨0.1⟩ s počáteční podmínkou y(0)=1

Definujme funkci f(t,y)=t y

f[t_,y_] :=t y ;

Řešme diferenciální rovnici  metodou Runge-Kutta pomocí našeho programu RK4. Zvolme ktok dt=0.2.

r=RK4[f,0,1,1.0,0.2,1]

Znázorníme si řešení graficky

g1=ListPlot[r,PlotStyle->PointSize[0.015]]

Nyní použijeme vestavěnou funkci pro řešení diferenciální rovnice.

res=NDSolve[{y'[t]==f[t,y[t]],y[0]==1},y,{t,0,1}]

Výsledek jsme dostali jako interpolační funkci z vypočtených dat. Data byla vypočtena metodou , kterou vybere program Mathematica jako vhodnou. Chceme-li libovolnou hodnotu řešení stačí napsat:

y[0.2]/.res

Chceme-li znázornit řešení graficky, stačí napsat následující příkaz:

g2=Plot[y[t]/.res,{t,0,1}]

Nyní porovnáme graficky oba výsledky:

Show[{g1,g2}]

Vypočteme přesné řešení:

rs=DSolve[{y'[t]==f[t,y[t]],y[0]==1},y[t],t]

Nyní porovnáme  oba výsledky se skutečným řešením normou rozdílu:

v1=Table[r[[i,2]],{i,1,6}];
v2=Table[y[(i-1) 0.2]/.res[[1]],{i,1,6}];
v=Table[Exp[((i-1) 0.2)^2/2.],{i,1,6}]

Chyba našeho výsledku

Norm[v1-v]/Sqrt[6]

Chyba  výsledku získaného programem NDSolve

Norm[v2-v]/Sqrt[6]

Vydíme, že jsme naším programem dostali pomerně dobrý výsledek i s krokem dt=0.2 .

Nyní zvolíme menší krok dt=0.01 a opět porovnáme  výsledek se skutečným výsledkem.

r2=RK4[f,0.0,1.0,1.0,0.01,0];

v1=Table[r2[[i,2]],{i,1,101}];
v2=Table[y[(i-1) 0.01]/.res[[1]],{i,1,101}];
v=Table[Exp[((i-1) 0.01)^2/2.],{i,1,101}];

Norm[v1-v]/Sqrt[101]

Norm[v2-v]/Sqrt[101]

Dostali jsme lepší výsledek než programem NDSolve. To je způsobeno tím, že program NDSolve nepožaduje větší přesnost než (10^(-7)). To ale můžeme ovlivnit pomocí options.

res=NDSolve[{y'[t]==f[t,y[t]],y[0]==1},y,{t,0,1},AccuracyGoal->14,PrecisionGoal->14,Method->ExplicitRungeKutta]

v1=Table[r2[[i,2]],{i,1,101}];
v2=Table[y[(i-1) 0.01]/.res[[1]],{i,1,101}];
v=Table[Exp[((i-1) 0.01)^2/2.],{i,1,101}];

Norm[v1-v]/Sqrt[101]

Norm[v2-v]/Sqrt[101]

Program RKSolve pro výpočet řešení  diferenciální rovnice y'=f(t,y) na intervalu ⟨(Subscript[t, 0]),(Subscript[t, 1])⟩ s počáteční podmínkou y((Subscript[t, 0]))=(Subscript[y, 0]) s krokem dt metodou Runge-Kutta. Zde jsou použity složitější programovací prvky. Tento program řeší i soustavy diferenciálních rovnic.

RKStep[f_, y_, y0_, dt_] :=
	Module[{ k1, k2, k3, k4 },
		k1 = dt N[ f /. Thread[y -> y0] ];
		k2 = dt N[ f /. Thread[y -> y0 + k1/2]];
		k3 = dt N[ f /. Thread[y -> y0 + k2/2]];
		k4 = dt N[ f /. Thread[y -> y0 + k3] ];
		y0 + (k1 + 2 k2 + 2 k3 + k4)/6
	]
RKSolve[f_List, y_List, y0_List, {t1_,n_}] :=
	NestList[ RKStep[f, y, #, N[t1/n]]&, N[y0],n ];
RKSolve[f_List, y_List, y0_List, {t_, t0_, t1_, n_}] :=
	Module[{res},
		res = RKSolve[ Prepend[f,1],Prepend[y,t],Prepend[y0,t0], {t1 - t0,n} ];
 res]

Mějme stejnou rovnici y'=t y na intervalu ⟨0,1⟩ s počáteční podmínkou y(0)=1. Zvolíme delení intervalu na 5 dílků, tedy dt=0.2.

f[t_,y_]=y  t;

data=RKSolve[{f[t,y]},{y},{1.0},{t,0,1,5}]

ListPlot[data,PlotStyle->PointSize[0.015]]

Příklad 2:
Mějme soustavu rovnic x'= y  y'= x y na intervalu ⟨0,1⟩ s počáteční podmínkou x(0)=1, y(0)=1. Zvolíme delení intervalu na 5 dílků.

f1[t_,x_,y_]=y;
f2[t_,x_,y_]=x y;

RKSolve[{f1[t,x,y],f2[t,x,y]},{x,y},{1.,1.},{t,0,1,5}]

Příklad 3:
Mějme soustavu rovnic x'= y  y'=- x na intervalu ⟨0,2π⟩ s počáteční podmínkou x(0)=1, y(0)=1. Zvolíme dělení intervalu na 100 dílků.

f1[t_,x_,y_]=y ;
f2[t_,x_,y_]=- x;

data=RKSolve[{f1[t,x,y],f2[t,x,y]},{x,y},{1.0,1.0},{t,0,2 Pi,100}];

Vykreslíme si křivku zadonou parametrickými rovnicemi
x=x(t)
y=y(t)
kde x(t) a y(t) jsou řešení soustavy rovnic.

data3=Drop[#,1]& /@ data;

ListPlot[data3,PlotStyle->PointSize[0.02],AspectRatio->1]

Nyní vykreslíme řešení x(t) v závislosti na t.

data1=Drop[#,-1]& /@ data;

ListPlot[data1,PlotStyle->PointSize[0.01]]

Nyní vykreslíme řešení y(t) v závislosti na t.

data2=Drop[#,{2,2}]& /@ data;

ListPlot[data2,PlotStyle->PointSize[0.01]]

Příklad 4:
Mějme soustavu rovnic x'= (x^2)+(y^2)-9  y'=(x^2)-1 na intervalu ⟨0,5⟩ s počáteční podmínkou x(0)=-0.5, y(0)=1.2 . Zvolíme dělení intervalu na 100 dílků. Vykreslíme příslušné grafy.

f1[t_,x_,y_]=x x+y y-9 ;
f2[t_,x_,y_]=x x-1;

data=RKSolve[{f1[t,x,y],f2[t,x,y]},{x,y},{-0.5,1.2},{t,0,5,100}];

data3=Drop[#,1]& /@ data;

ListPlot[data3,PlotStyle->PointSize[0.015]]

data1=Drop[#,-1]& /@ data;

ListPlot[data1,PlotStyle->PointSize[0.017],PlotRange->All]

data2=Drop[#,{2,2}]& /@ data;

ListPlot[data2,PlotStyle->PointSize[0.017],PlotRange->All]

Program RK4Plus pro výpočet řešení  diferenciální rovnice y'=f(t,y) na intervalu ⟨(Subscript[t, 0]),(Subscript[t, 1])⟩ s počáteční podmínkou y((Subscript[t, 0]))=(Subscript[y, 0]) s automatickou volbou kroku metodou Runge-Kutta. Tento program řeší i soustavy diferenciálních rovnic.

RKStep[f_, y_, y0_, dt_] :=
	Module[{ k1, k2, k3, k4 },
		k1 = dt N[ f /. Thread[y -> y0] ];
		k2 = dt N[ f /. Thread[y -> y0 + k1/2]];
		k3 = dt N[ f /. Thread[y -> y0 + k2/2]];
		k4 = dt N[ f /. Thread[y -> y0 + k3]];
		y0 + (k1 + 2 k2 + 2 k3 + k4)/6
	];
RK4Plus[f_,y_,t_,t0_,t1_, y0_, dt_,dtmin_,accur_] :=
	Module[{ta,ya,h1,ras1,res2,respruh,resplus,s ,m,vys},
		ta=t0;
		ya=y0;
		h=dt;
		m=0;
		vys={Prepend[ya,ta]};
	  	While[Abs[t1-ta]>0.0000001 && m==0,
	     	If[Abs[t1-ta]<Abs[h],h=t1-ta];
			   res1=RKStep[Prepend[f,1],Prepend[y,t],Prepend[ya,ta],h];
		     s=999999;
			  While[s>accur && m==0,
			  respruh=RKStep[Prepend[f,1],Prepend[y,t],Prepend[ya,ta],h/2];
			  res2=RKStep[Prepend[f,1],Prepend[y,t],Prepend[Drop[respruh,1],ta+h/2],h/2];
			  resplus=(16*res2-res1)/15;
				s=Sqrt[(Drop[res2,1]-Drop[resplus,1]).(Drop[res2,1]-Drop[resplus,1])];
				If[s>accur && m==0,
					  If[h>dtmin,
				         Print[ "t=",ta,"  Pulim krok  z ",
							         "h=",h, "  s=",s];h=h/2;res1=respruh,
					      Print["t=",ta+h,"   h=",h," je mensi nez hmin   ", "t=",ta,"  s=",s];m=1;ta=ta+h],
					 vys=Join[vys,{resplus}];
						If[s<0.8*accur/32,
						       Print["t= ",ta+h,"   h=",h,"  Zdvojnasobuji krok    ","  s=",s];h1=2*h;m=0,h1=h];
		         	    ya=Drop[resplus,1];ta=ta+h;h=h1];
			];
		];
 N[vys,15]
	]

Příklad 5:
Mějme  rovnici   y'= sin(3π x y) na intervalu ⟨0.1⟩ s počáteční podmínkou x(0)=1, y(0)=1. Zvolíme automatickou volbu kroku..

f[x_,y_]=Sin[3 Pi y   x];

data=RK4Plus[{f[x,y]},{y},x,0.0,3.0,{1.0},0.1,0.001,0.000001]

ListPlot[data,PlotStyle->PointSize[0.02]]