Способы решения линейных дифференциальных уравнений (ЛДУ), страница 2

 x₃^’(t)= y^’’’(t)=(10-5* y^’-106* y^’’- y)/12=(10-5*x₂-106*x₃-x₁)/12

3. Текст программы:

mf.m

//Вычисляет правую часть дифференциального уравнения для шага t, решения y

function z=mf(t,y,a,f)

z=a*y+f;

oscil.m

//Задает уравнение для солвера

function f=oscil(t,y)

f=[y(2);y(3);(10-y(1)-106*y(3)-5*y(2))/12];

RK.m

//Формула Рунге-Кутта 4-го порядка для расчета начальных точек

function y=RK(t,y,a,f,h)

    k1=mf(t,y,a,f);

    k2=mf(t+h/2,y+h/2*k1,a,f);

    k3=mf(t+h/2,y+h/2*k2,a,f);

    k4=mf(t+h,y+h*k3,a,f);

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

TR.m

//Главная функция, осуществляющая аналитическое решение, решение при помощи солвера ode45 и численный метод Адамса-Башфорта 3-го порядка.

%***********************%

%Exact Analytic Solution%

%***********************%

syms y

digits(20);

y =vpa( dsolve('12*D3y+106*D2y+5*Dy+y = 10', 'y(0)=0', 'D2y(0) = 0','Dy(0)=0', 'x'));

figure;

ezplot(y,[0,200])

s=vectorize(y)

x=0:200;

title('Exact Analytic Solution')

%**********%

%A-B method%

%**********%

%Define constants

h=0.05; DT=200; n=3;

a=[0 1 0;0 0 1;-1/12 -5/12 -106/12];

f=[0;0;10/12];

y=[0;0;0];t=0;

%end of definition

figure;

m=floor(DT/h);

out=zeros(m,n+1);

%Set 3 first point for starting A-B method, using R-K method

%%FIRST%%

p=[t,y'];

out(1,:)=p;%

k1=y;

f1=mf(t,k1,a,f);

t=t+h;

%%SECOND%%

k2=RK(t,k1,a,f,h);

p=[t,k2'];

out(2,:)=p;

f2=mf(t,k2,a,f);

t=t+h

%%THIRD%%

k3=RK(t,k2,a,f,h);

p=[t,k3'];

out(3,:)=p;

f3=mf(t,k3,a,f);

t=t+h;

for i=4:m; %A-B method start

    y=k3+h*(23*f3-16*f2+5*f1)/12;