У табл. 2.1 наведено приклади програм для наближеного розв’язання диференціальних рівнянь за допомогою вище розглянутих однокрокових і багатокрокових числових методів, написаних мовою пакета Matlab.

Однокрокові методи

Метод Ейлера

Метод Гюна

function E=euler(a,b,ya,h)

M=(b-a)/h;

t=zeros(1,M+1);

y=zeros(1,M+1);

t=a:h:b;

y(1)=ya;

for j=1:M

y(j+1)=y(j)+h*feval('f',t(j),y(j));

end

E=[t'y'];

function f=f(t,y)

f=0.5*(t-y);

function H=heun(a,b,ya,h)

M=(b-a)/h;

t=zeros(1,M+1);

y=zeros(1,M+1);

t=a:h:b;

y(1)=ya;

for j=1:M

k1=feval('f',t(j),y(j));

k2=feval('f',t(j+1),y(j)+h*k1);

y(j+1)=y(j)+h*feval('f',t(j),y(j));

end

H=[t'y'];

function f=f(t,y)

f=0.5*(t-y);

Метод рядів Тейлора 4-го порядку

Метод Рунге-Кутта 4-го порядку

function T4=taylor(a,b,ya,h)

M=(b-a)/h;

t=zeros(1,M+1);

y=zeros(1,M+1);

t=a:h:b;

y(1)=ya;

for j=1:M

D=feval('df',t(j),y(j));

y(j+1)=y(j)+h*(D(1)+h*(D(2)/2+h*

(D(3)/6+h*D(4)/24)));

end

T4=[t'y'];

function df=df(t,y)

df=[(t-y)/2 (2-t+y)/4 (-2+t-y)/8 (2-t+y)/16];

function R=rk4(a,b,ya,h)

M=(b-a)/h;

t=zeros(1,M+1);

y=zeros(1,M+1);

t=a:h:b;

y(1)=ya;

for j=1:M

k1=h*feval(‘f’,t(j),y(j));

k2=h*feval(‘f’,t(j)+h/2,y(j)+k1/2);

k3=h*feval(‘f’,t(j)+h/2,y(j)+k2/2);

k4=h*feval(‘f’,t(j)+h,y(j)+k3);

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

end

R=[t'y'];

function f=f(t,y)

f=0.5*(t-y);

Багатокрокові методи

Метод Адамса-Бешфорса-Маултона

Метод Мілна-Сімпсона

function A=abm(a,b,ya,h)

M=(b-a)/h;

t=zeros(1,M+1);

y=zeros(1,M+1);

t=a:h:b;

n=length(t);

R=rk4(a,b,ya,h);

for j=1:4

y(1,j)=R(j,2);

end

if n<5,break,end;

F=zeros(1,4);

function M=milne(a,b,ya,h)

M=(b-a)/h;

t=zeros(1,M+1);

y=zeros(1,M+1);

t=a:h:b;

n=length(t);

R=rk4(a,b,ya,h);

for j=1:4

y(1,j)=R(j,2);

end

if n<5,break,end;

F=zeros(1,4);

F=feval('f',t(1:4),y(1:4));

for k=4:n-1

p=y(k)+(h/24)*(F*[-9 37 -59 55]');

t(k+1)=t(1)+h*k;

F=[F(2) F(3) F(4) feval('f',t(k+1),p)];

y(k+1)=y(k)+(h/24)*(F*[1 -5 19 9]');

F(4)=feval('f',t(k+1),y(k+1));

end

A=[t'y'];

function f=f(t,y)

f=0.5*(t-y);

F=feval('f',t(1:4),y(1:4));

pold=0;

yold=0;

for k=4:n-1

pnew=y(k-3)+(4*h/3)*(F(2:4)*[2 -1 2]');

pmod=pnew+28*(yold-pold)/29;

t(k+1)=t(1)+h*k;

F=[F(2) F(3) F(4) feval('f',t(k+1),pmod)];

y(k+1)=y(k-1)+(h/3)*(F(2:4)*[1 4 1]');

pold=pnew;

yold=y(k+1);

F(4)=feval('f',t(k+1),y(k+1));

end

M=[t'y'];

function f=f(t,y)

f=0.5*(t-y);

Метод Хемінга

function H=hamming(a,b,ya,h)

M=(b-a)/h;

t=zeros(1,M+1);

y=zeros(1,M+1);

t=a:h:b;

n=length(t);

R=rk4(a,b,ya,h);

for j=1:4

y(1,j)=R(j,2);

end

if n<5,break,end;

F=zeros(1,4);

F=feval('f',t(1:4),y(1:4));

pold=0;

cold=0;

for k=4:n-1

pnew=y(k-3)+(4*h/3)*(F(2:4)*[2 -1 2]');

pmod=pnew+112*(cold-pold)/121;

t(k+1)=t(1)+h*k;

F=[F(2) F(3) F(4) feval('f',t(k+1),pmod)];

cnew=(9*y(k)-y(k-2)+3*h*(F(2:4)*[-1 2 1]'))/8;

y(k+1)=cnew+9*(pnew-cnew)/121;

pold=pnew;

cold=cnew;

F(4)=feval('f',t(k+1),y(k+1));

end

H=[t'y'];

function f=f(t,y)

f=0.5*(t-y);