Методи

explicitEuler[f_,h_,t0_,T_,u0_]:= (
n=Ceiling[(T-t0)/h];
t=Table[t0+i*h,{i,0,n}];
y=Table[0,{n+1}];

y[[1]] =u0;
For[i=1, i<n+1, i++,
y[[i+1]] = y[[i]]+h*f[t[[i]],y[[i]]]
];

Transpose[{t,y}]
)

implicitEuler[f_,h_,t0_,T_,u0_]:= (
n=Ceiling[(T-t0)/h];
t=Table[t0+i*h,{i,0,n}];
y=Table[0,{n+1}];

y[[1]] =u0;
For[i=1, i<n+1, i++,
initial = y[[i]]+h*f[t[[i]],y[[i]]];
y[[i+1]]=yNew/.FindRoot[yNew== y[[i]]+h*f[t[[i+1]],yNew],{yNew,initial}];
];
Transpose[{t,y}]
)

RK4[f_,h_,t0_,T_,u0_]:= (
n=Ceiling[(T-t0)/h];
t=Table[t0+i*h,{i,0,n}];
y=Table[0,{i, 0, n}];
y[[1]] = u0;
For[i=1, i<=n, i++,
k1 = h*f[t[[i]],y[[i]]];
k2 = h*f[t[[i]]+h/2,y[[i]]+k1/2];
k3 = h*f[t[[i]]+h/2,y[[i]]+k2/2];
k4 =h* f[t[[i]]+h,y[[i]]+k3];
y[[i+1]] = y[[i]]+1/6k1+1/3k2+1/3k3+1/6k4];

Transpose[{t,y}]
)

AB4[f_,h_,t0_,T_,u0_]:=(
n=Ceiling[(T-t0)/h];
t=Table[t0+i h,{i,0,n}];
y=Table[0,{n+1}];
fi=Table[0,{n+1}]; (*list of values of the right-hand side f*)
y[[1]]=u0;
(*Compute the first approximate values, using the RK4 method*)
For[i=1, i≤3, i++,
k1=h f[t[[i]],y[[i]]];
k2=h f[t[[i]]+h/2,y[[i]]+k1/2];
k3=h f[t[[i]]+h/2,y[[i]]+k2/2];
k4=h f[t[[i]]+h,y[[i]]+k3];
y[[i+1]]=y[[i]]+(1/6)(k1+2k2+2k3+k4);
fi[[i+1]]=f[t[[i+1]],y[[i+1]]];

];
(*Compute the remaining values, using the AB4 formula*)
For[i=4,i<n+1,i++,
y[[i+1]]=y[[i]]+((h )/24)(55 fi[[i]]- 59 fi[[i-1]]+ 37 fi[[i-2]]- 9fi[[i-3]]);
fi[[i+1]]=f[t[[i+1]],y[[i+1]]];
];
Transpose[{t,y}]
)

НЕ Е РЕАЛИЗИРАНО НА УПРАЖНЕНИЕ!
Тяло, падащо свободно към земната повърхност, може да бъде моделирано чрез системата ОДУ от втори ред:
x''(t)=9.8, 0<t<=10,
x(0) = 0 //началното положение е 0
x’(0) = 0 //началната скорост, с която е пуснато тялото е 0.

Свеждаме задачата към задачата на Коши с полагането x'(t)=v. Тогава:
x'=v
v'=9.8
x(0)=0
v(0)=0
или 
Свеждаме задачата към задачата на Коши с полагането x'(t)=v. Тогава:
u'=f(t,u)
u(0)=(Subscript[u, 0])
където u=(x,v)

Clear[x,t];
exact1[t_]=x[t]/.DSolve[{x''[t]==9.8,x[0]==0, x'[0]==0},x[t],t][[1]];
h=0.01;
t=Table[0+i*h,{i,0,10/h}];

(*Ще модифицираме кода за метода на Ойлер да работи за 2 уравнения*)
explicitEuler[f1_,f2_,h_,t0_,T_,x0_,v0_]:= (
n=Ceiling[(T-t0)/h];
t=Table[t0+i*h,{i,0,n}];
x=Table[0,{n+1}];
v=Table[0,{n+1}];

x[[1]] =x0;
v[[1]]=v0;
For[i=1, i<n+1, i++,
x[[i+1]] = x[[i]]+h*f1[t[[i]],x[[i]],v[[i]]]; 
v[[i+1]] = v[[i]]+h*f2[t[[i]],x[[i]],v[[i]]]
];

Transpose[{t,x,v}]
)

f1[t_,x_,v_]:=v;
f2[t_,x_,v_]:=9.8;
apprPos1 = explicitEuler[f1,f2,h,0,10,0,0][[All, 2]];
Animate[Graphics[{Circle[{0,-exact1[t[[i]]]},5], {Red,Dashed,Circle[{0,-apprPos1[[i]]},5]}},
							 PlotRange->{{0,10},{-550,0}}, Axes->{False,True}],
		 {i,1,Length[t],1}]

apprPos1[[90]]

explicitEuler[f1_,f2_,h_,t0_,T_,x0_,v0_]:= (
n=Ceiling[(T-t0)/h];
t=Table[t0+i*h,{i,0,n}];
x=Table[0,{n+1}];
v=Table[0,{n+1}];

x[[1]] =x0;
v[[1]]=v0;
For[i=1, i<n+1, i++,
x[[i+1]] = x[[i]]+h*f1[t[[i]],x[[i]],v[[i]]]; 
v[[i+1]] = v[[i]]+h*f2[t[[i]],x[[i]],v[[i]]]
];

Transpose[{t,x,v}]
)

(*Това модифициране на кода по-горе не е извънредно удобно, защото налага да пренаписваме кода за различен брой уравнения. Освен това реално кодът маркиран в жълто прави едно и също нещо за x и v. Тогава реално y може да е списък от вида {{_,_},{_,_},...,{_,_}}. Тогава f(t, u)=(((f1(t,u), f2(t,u)))^T) *)
f[x_,u_]:={u[[2]],9.8} 
explicitEuler[f_,h_,t0_,T_,u0_]:= (
n=Ceiling[(T-t0)/h];
t=Table[t0+i*h,{i,0,n}];
y=Table[0,{n+1}];

y[[1]] =u0;
For[i=1, i<n+1, i++,
y[[i+1]] = y[[i]]+h*f[t[[i]],y[[i]]]
];

Transpose[{t,y}]
)

apprPos2= explicitEuler[f,h,0,10,{0,0}][[All, 2, 1]];
apprPos2-apprPos1 (*получените резултати с двата метода са едни и същи*)

Animate[Graphics[{Circle[{0,-exact1[t[[i]]]},5], {Red,Dashed,Circle[{0,-apprPos2[[i]]},5]}},
							 PlotRange->{{0,10},{-550,0}}, Axes->{False,True}],
		 {i,1,Length[t],1}]

Решете задачата
r'(s)=cos(phi), 
z'(s)=sin(phi), 
phi'(s)=2b+ c z -((sin(phi))/r), 
r(0) =z(0)=phi(0)= 0,
използвайки Адамс-Башфорт. Решавайте, докато r<1

b=1/0.5424;
c=-2.9;

f[s_,u_]:={
Cos[u[[3]]],
Sin[u[[3]]],
If[u[[1]]≠0,2b+c*u[[2]]-Sin[u[[3]]]/u[[1]],b]
}
res=AB4[f,0.01,0,5,{0,0,0}];
ListPlot[Table[{res[[i,2,1]],res[[i,2,2]]},{i,1,Length[res]}], AspectRatio->Automatic]

(*Мерните единици в задачата са избрани обикновено така, че радиусът на капиляра е 1. Това означава, че трябва да решаваме, докато r(s)≤1. След това резултатите, които се получават, не са физически смислени. За целта модифицираме функцията AB4*)
AB4[f_,h_,t0_,u0_]:=(
t={t0};
y={u0};
fi={f[t0,u0]}; (*list of values of the right-hand side f*)
(*Compute the first approximate values, using the RK4 method*)
For[i=1, i≤3, i++,
AppendTo[t,t0+i*h];
k1=h f[t[[i]],y[[i]]];
k2=h f[t[[i]]+h/2,y[[i]]+k1/2];
k3=h f[t[[i]]+h/2,y[[i]]+k2/2];
k4=h f[t[[i]]+h,y[[i]]+k3];
AppendTo[y,y[[i]]+(1/6)(k1+2k2+2k3+k4)];
AppendTo[fi, f[t[[i+1]],y[[i+1]]]];
];

(*Compute the remaining values, using the AB4 formula*)
i=4;
While[y[[i,1]]≤1,
AppendTo[t,t0+i*h];
AppendTo[y,y[[i]]+((h )/24)(55 fi[[i]]- 59 fi[[i-1]]+ 37 fi[[i-2]]- 9fi[[i-3]])];
AppendTo[fi,f[t[[i+1]],y[[i+1]]]];
i++;
];
Transpose[{t,y}]
)

b=1/0.5424;
c=-2.9;

f[s_,u_]:={
Cos[u[[3]]],
Sin[u[[3]]],
If[u[[1]]≠0,2b+c*u[[2]]-Sin[u[[3]]]/u[[1]],b]
}
res=AB4[f,0.01,0,{0,0,0}];
ListPlot[Table[{res[[i,2,1]],res[[i,2,2]]},{i,1,Length[res]}], AspectRatio->Automatic]

Решете задачата
u'(x)=(10^6)(-u+sin t ), 0<t<=10,
u(0) = 0,
използвайки явен и неявен метод на Ойлер, Рунге-Кута и Адамс-Башфорт. Сравнете резултатите.

Clear[x];
f[x_,u_]:=10^6(-u+Sin[x])
exact3[x_,l_]:=l*Exp[-l*x]/(1+l^2)+l^2/(1+l^2)Sin[x]-l/(1+l^2)Cos[x];
plots[appr_]:=If[Max[Abs[appr]]<10^6,
				Show[
					ListPlot[appr,PlotRange->Automatic],
					Plot[exact3[x,10.^6],{x,0,10}, PlotStyle->Red]
				],
				"The step-size is too big!"]

h=0.1; (*За неявния метод на Ойлер получаваме добри резултати дори за голяма стъпка h=0.1. Причината за това, че методът е абсолютно устойчив, т.е. можем да разчитаме, че дори при големи стъпки няма да има лошо поведение на метода. *)
yImplicit = implicitEuler[f,h,0.,10.,0.];
plots[yImplicit]

(*Всички явни методи се провалят, защото те са условно устойчиви. За да получим добри резултатите стъпката трябва да е много малка, което забавя значително изчисленията.*)

h=0.00001;
yExplicit = explicitEuler[f,h,0.,10.,0.];
plots[yExplicit]

h=0.001;
yRK4 = RK4[f,h,0.,10.,0.];
plots[yRK4]

h=0.00001;
yAB4 = AB4[f,h,0.,10.,0.];
plots[yAB4]