Решете моделната задача
u’(x)=-10u, 0<x<=1,
u(0)=1
при n=4, 20, 100 и сравнете с точното решение.

Явен метод на Ойлер

explicitEuler[n_,f_]:= (
h=1.0/n;(*compute numerically*)
xList=Table[i*h,{i,0,n}];
yList=Table[0,{n+1}];

yList[[1]] = 1;
For[i=1, i<n+1, i++,
yList[[i+1]] = yList[[i]]+h*f[xList[[i]],yList[[i]]]
];

Transpose[{xList,yList}]
)

(*Compute the exact solution*)
f[x_,u_]:=-10*u
exact[t_]=DSolve[{u'[t]==f[t,u[t]],u[0]==1},u[t],t][[1,1,2]];
(*exact[x_]:=E^(-10x);*)
plotExact = Plot[exact[x],{x,0, 1}, PlotStyle->Red, PlotRange->All];

n=100(*4,20*);
appr = explicitEuler[n,f];
plotAppr = ListPlot[appr, PlotMarkers->●];
Show[plotExact, plotAppr]

Manipulate[Show[ListLinePlot[explicitEuler[n,f], PlotMarkers->●, PlotRange->{{0,1},{-2,2}}], plotExact],{n, 4, 100,4}]

Неявен метод на Ойлер

implicitEuler[n_,f_]:= (
h=1.0/n;
xList=Table[i*h,{i,0,n}];
yList=Table[0,{n+1}];

yList[[1]] = 1;
For[i=1, i<n+1, i++,
initialGuess =yList[[i]]+h*f[xList[[i]],yList[[i]]];
yList[[i+1]]=FindRoot[yNew==  yList[[i]]+h*f[xList[[i+1]],yNew],{yNew,initialGuess}][[1,2]]
];

Transpose[{xList,yList}]
)

(*Compute the exact solution*)
f[x_,u_]:=-10*u
exact[t_]=DSolve[{u'[t]==f[t,u[t]],u[0]==1},u[t],t][[1,1,2]];
plotExact = Plot[exact[x],{x,0, 1}, PlotStyle->Red, PlotRange->All];

Manipulate[Show[ListLinePlot[implicitEuler[n,f], PlotMarkers->●, PlotRange->{{0,1},{-2,2}}], plotExact],{n, 4, 100,4}]

Неявен метод на Ойлер – решаване на нелинейните уравнения с метод на Нютон

(*Newton method solves equations g(x)=0, generating a sequence of points 
x0,x1,x2,...,xn,....
using the formula
(Subscript[x, n])=(Subscript[x, n-1])-((g[(Subscript[x, n-1])))/(g'[(Subscript[x, n-1]))))
The method stops if 
the difference between two consequtive iterations is less than tol
the number of iterations is greater than the maxIter
*)

Newton[g_,x0_, tol_, maxIter_]:=(
xCurrent = x0;
xNext = xCurrent - g[xCurrent]/g'[xCurrent];
iter=1;

While[Abs[xCurrent-xNext]≥tol && iter<maxIter,
xCurrent = xNext;
xNext = xCurrent - g[xCurrent]/g'[xCurrent];
iter++;
];

xNext
)

g[x_]:=x^2-2x+1.+Sin[x]+Log[x]
Newton[g,10,0.0000001,100]

FindRoot[g[x]==0,{x,10}]

implicitEulerWithNewton[n_,f_]:= (
h=1.0/n;
xList=Table[i*h,{i,0,n}];
yList=Table[0,{n+1}];

yList[[1]] = 1;
For[i=1, i<n+1, i++,
initialGuess =yList[[i]]+h*f[xList[[i]],yList[[i]]];
g[yNew_]:=yNew-(  yList[[i]]+h*f[xList[[i+1]],yNew]);
yList[[i+1]]=Newton[g,initialGuess,0.0000001,100];
];

Transpose[{xList,yList}]
)

f[x_,u_]:=-10*u
exact[t_]=DSolve[{u'[t]==f[t,u[t]],u[0]==1},u[t],t][[1,1,2]];
plotExact = Plot[exact[x],{x,0, 1}, PlotStyle->Red, PlotRange->All];

Manipulate[Show[ListLinePlot[implicitEulerWithNewton[n,f], PlotMarkers->●, PlotRange->{{0,1},{-2,2}}], plotExact],{n, 4, 100,4}]