Method Runge-Kutta

ClearAll["Global`*"];

Comparing the Runge-Kutta method with the NDSolve method

Another part of the CAO-RungeKutta.nb notebook shows the comparison of the results of the derived method with NDSolve. We recommend to play with different values of the number of steps n, or change the specified function f (myf is called myf). Note, in particular, what the jump problem is in a derivative (why myf is entered with If).Where is a discontinuity f, there are unreliable numerical solutions of differential equations. Not only there ... fortunately, These function are not in the electrical circuit theory, because nature forbids these type of function.

tStart=0.;
tEnd=15;
yStart=5;
myf[t_,y_]:=If[t<0.5*(tStart+tEnd),Sin[3t-.5y],-y*Cos[t+(y^2)]];
n=250;
Δt=((tEnd-tStart)/n);

res=NestList[dal,{tStart,yStart},n];
resNDS=NDSolve[{y'[t]==myf[t,y[t]],y[tStart]==yStart},y,{t,tStart,tEnd}][[1]];

Using the Runge-Kutta methods 4th order

... In the next section it is shown using the Runge-Kutta methods fourth order (coincided with the five members of the Taylor series, i.e. up to the fourth power including t). ...

derivace[t_,y_]:=myf[t,y];
RungeKuttaKrok[{t_,y_}]:=Module[
{k1,k2,k3,k4},
k1=derivace[t,y];
k2=derivace[t+((Δt)/2),y+((Δt)/2)*k1];
k3=derivace[t+((Δt)/2),y+((Δt)/2)*k2];
k4=derivace[t+Δt,y+Δt*k3];
{t+Δt,y+((Δt)/6)*(k1+2k2+2k3+k4)}
];
res2=NestList[RungeKuttaKrok,{tStart,yStart},n];
Show[Plot[y[t]/.resNDS,{t,tStart,tEnd},PlotRange->All,PlotStyle->Thickness[0.005],AxesLabel->{"t[s]","y"}],ListPlot[res2,PlotStyle->{Red,PointSize[0.01]}]]

Graphic demonstration methods of Runge-Kutta 4th order

We have to solve differential equation y’[t]==-Sin[(t^3)+y[t]+0.8]

der[{t_,y_}]:=-Sin[(t^3)+y+0.8];

tStart=1.5;
yStart=0.5;
Δt=0.35;

First, we show the solution “exactly” calculated by the NDSolve method (calculated in any case with a step much smaller than  Δt):

rNDS=NDSolve[{y'[t]==der[{t,y[t]}],y[tStart]==yStart},y,{t,tStart,tStart+Δt}][[1]];
plNDS=Plot[y[t]/.rNDS,{t,tStart,tStart+Δt},Ticks->{{{tStart,"t"},{tStart+((Δt)/2),"t+((Δt)/2)"},{tStart+Δt,"t+Δt"}}},AxesOrigin->{tStart-0.2*Δt,0.4},PlotRange->{{tStart-0.2*Δt,tStart+1.2*Δt},{0.4,0.8}},AxesLabel->{"t","y"},PlotStyle->{Thickness[0.0055],Black}]

We evaluate the k1, derivative at the beginning (red point), i.e. at [t, y (t)]. The derivative k1 isthe slope of the red dashed line.

bodStart={tStart,yStart};
znackaBodStart=Graphics[{Red,PointSize[0.015],Point[bodStart]}];
k1=der[bodStart];
plk1=Plot[yStart+(t-tStart)*k1,{t,tStart-0.2Δt,tStart+Δt},PlotStyle->{Red,Dashed,Thickness[0.004]}];
Show[plNDS,znackaBodStart,plk1]

We will start from starting point with sope k1 (red arrow).
We use step ((Δt)/2) and reach the point [t+((Δt)/2), y(t) + k1 ((Δt)/2)] (blue point).

bodProk2={tStart+((Δt)/2),yStart+k1*((Δt)/2)};
plk1krok=Graphics[{Red,Arrow[{bodStart,bodProk2}]}];
znackaBodProk2=Graphics[{Blue,PointSize[0.015],Point[bodProk2]}];
Show[plNDS,znackaBodStart,plk1,plk1krok,znackaBodProk2]

We evaluate value of the derivative at the blue point k2. Derivative k2 is the slope of dashed blue line.

k2=der[bodProk2];
plk2=Plot[bodProk2[[2]]+(t-tStart-((Δt)/2))*k2,{t,tStart,tStart+Δt},PlotStyle->{Blue,Dashed,Thickness[0.004]}];
Show[plNDS,znackaBodStart,plk1,plk1krok,znackaBodProk2,plk2]

We start again from the red point, but in the direction with slope  k2, that is parallel with the blue line (blue arrow). 
The step is ((Δt)/2), and we reach  the point [t+((Δt)/2), y(t) + k2 ((Δt)/2)], green point).

bodProk3={tStart+((Δt)/2),yStart+k2*((Δt)/2)};
plk2krok=Graphics[{Blue,Arrow[{bodStart,bodProk3}]}];
znackaBodProk3=Graphics[{Green,PointSize[0.015],Point[bodProk3]}];
Show[plNDS,znackaBodStart,plk1,plk1krok,znackaBodProk2,plk2,plk2krok,znackaBodProk3]

We evaluate new value of derivative at the green point k3.  k3 is the slope of green dashed line k3.

k3=der[bodProk3];
plk3=Plot[bodProk3[[2]]+(t-tStart-((Δt)/2))*k3,{t,tStart,tStart+Δt},PlotStyle->{Green,Dashed,Thickness[0.004]}];
Show[plNDS,znackaBodStart,plk1,plk1krok,znackaBodProk2,plk2,plk2krok,znackaBodProk3,plk3]

Using  derivative  k3 and step  Δt we reach  point [t+((Δt)/2), y(t) + k3 Δt]  (orange point).

bodProk4={tStart+Δt,yStart+k3*Δt};
plk3krok=Graphics[{Green,Arrow[{bodStart,bodProk4}]}];
znackaBodProk4=Graphics[{Orange,PointSize[0.015],Point[bodProk4]}];
Show[plNDS,znackaBodStart,plk1,plk1krok,znackaBodProk2,plk2,plk2krok,znackaBodProk3,plk3,plk3krok,znackaBodProk4]

k4 is the derivative at the orange point and it is shown as dashed orange line.

k4=der[bodProk4];
plk4=Plot[bodProk4[[2]]+(t-tStart-Δt)*k4,{t,tStart,tStart+1.2Δt},PlotStyle->{Orange,Dashed,Thickness[0.004]}];
Show[plNDS,znackaBodStart,plk1,plk1krok,znackaBodProk2,plk2,plk2krok,znackaBodProk3,plk3,plk3krok,znackaBodProk4,plk4]

We evaluate final derivative kVys as a weighted average from derivatives k1, k2, k3 a k4.

kVys=(1/6)*(k1+2k2+2k3+k4);

Using this derivative  kVys and Δt we reach the black point.

bodVys={tStart+Δt,yStart+kVys*Δt};
znackaBodVys=Graphics[{Black,PointSize[0.015],Point[bodVys]}];
plkVyskrok=Graphics[{Black,Arrow[{bodStart,bodVys}]}];
Show[plNDS,znackaBodStart,plk1,plk1krok,znackaBodProk2,plk2,plk2krok,znackaBodProk3,plk3,plk3krok,znackaBodProk4,plk4,plkVyskrok,znackaBodVys]

Nevertheles Δt is quite large, next value at time  t+Δt is relatively exact.

grafy={plNDS,znackaBodStart,plk1,plk1krok,znackaBodProk2,plk2,plk2krok,znackaBodProk3,plk3,plk3krok,znackaBodProk4,plk4,plkVyskrok,znackaBodVys};

kroky=Table[Take[grafy,i],{i,1,14}];

Manipulate[Show[kroky[[krok]]],{krok,Range[14],Setter}]