Методи

explicitEuler[f_,h_,t0_,T_,u0_]:= (
n=Ceiling[(T-t0)/h];
x=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[x[[i]],y[[i]]]
];

Transpose[{x,y}]
)

implicitEuler[f_,h_,t0_,T_,u0_]:= (
n=Ceiling[(T-t0)/h];
x=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[x[[i]],y[[i]]];
y[[i+1]]=yNew/.FindRoot[yNew== y[[i]]+h*f[x[[i+1]],yNew],{yNew,initial}];
];
Transpose[{x,y}]
)

RK4[f_,h_,t0_,T_,u0_]:= (
n=Ceiling[(T-t0)/h];
x=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[x[[i]],y[[i]]];
k2 = h*f[x[[i]]+h/2,y[[i]]+k1/2];
k3 = h*f[x[[i]]+h/2,y[[i]]+k2/2];
k4 =h* f[x[[i]]+h,y[[i]]+k3];
y[[i+1]] = y[[i]]+1/6k1+1/3k2+1/3k3+1/6k4];

Transpose[{x,y}]
)

AB4[f_,h_,t0_,T_,u0_]:=(
n=Ceiling[(T-t0)/h];
x=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[x[[i]],y[[i]]];
k2=h f[x[[i]]+h/2,y[[i]]+k1/2];
k3=h f[x[[i]]+h/2,y[[i]]+k2/2];
k4=h f[x[[i]]+h,y[[i]]+k3];
y[[i+1]]=y[[i]]+(1/6)(k1+2k2+2k3+k4);
fi[[i+1]]=f[x[[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[x[[i+1]],y[[i+1]]];
];
Transpose[{x,y}]
)