2020-10-10
159
произвольную точку определим
Система дифференциальных уравнений 
где x — независимая переменная, а y 1(x), y 2(x),..., yn (x) — неизвестные функции, n — порядок системы.
Обозначив 
или 
Решением системы называется вектор-функция 
, которая определена и непрерывно дифференцируема на интервале (a, b) и удовлетворяет системе, т.е. для всех 
справедливо 
Дифференциальные уравнения n-ого порядка 
Приводим к системе дифференциальных уравнений 
Программа. Дифференциальное уравнение. Метод Эйлера.
function DATA
global xn n x0 y0;
xn=3;
n=140;
x0=1;
y0=-2;
end
function [ dfdx ] = f(x,y)
dfdx=-y/x;
end
function [ y ] = Eiler(x,y,h)
y=y+h*f(x,y);
end
function [ x,y ] = fun_Dif_Ur_Eiler(x0,y0,xn,n)
h=(xn-x0)/n;
for i=1:n+1
if i==1
x(i)=x0;
y(i)=y0;
else
x(i)=x(i-1)+h;
y(i)=Eiler(x(i-1),y(i-1),h);
end %if
end %for i
end % function
function GLAV_Dif_Ur_Eiler
global xn n x0 y0 x y;
DATA;
[ x,y ] = fun_Dif_Ur_Eiler(x0,y0,xn,n);
REPORT;
end % function
function REPORT
global xn n x0 y0 x y;
disp('Diffferencialnoye uravneniye Euler');
disp(' ');
disp('Arguments');
disp(['x0 = ' num2str(x0,'%10.5f') ]);
disp(['y0 = ' num2str(y0,'%10.5f') ]);
disp(['xn = ' num2str(xn,'%10.5f') ]);
disp(['n = ' num2str(n,'%10.5f') ]);
disp('Results');
disp(' i x fx ');
|
|
|
disp(' ______________________________')
i=0;
for i=1:length(x)
xx=x(i);
ffx=y(i);
disp(sprintf('%10.3f\t%10.3f\t %10.3f',i,xx,ffx));
end %for
plot(x,y,'r.');
grid on;
xlabel('x');
ylabel('y');
title('Diffferencialnoye uravneniye Euler');
end % function
|
|
|
|
|
|
|
|
Программа. Дифференциальное уравнение.
Модифицированный метод Эйлера
function DATA
global xn n x0 y0;
xn=3;
n=14;
x0=1;
y0=-2;
end
function [ dfdx ] = f(x,y)
dfdx=-y/x;
end
function [ y ] = Modif_Eiler(x,y,h)
xnew=x+h/2;
ynew=y+h*f(x,y)/2;
y=y+h*f(xnew,ynew);
end
function [ x,y ] = fun_Dif_Ur_Modif_Eiler(x0,y0,xn,n)
h=(xn-x0)/n;
for i=1:n+1
if i==1
x(i)=x0;
y(i)=y0;
else
x(i)=x(i-1)+h;
y(i)=Modif_Eiler(x(i-1),y(i-1),h);
end %if
end %for i
end % function
function GLAV_Dif_Ur_Modif_Eiler
global xn n x0 y0 x y;
DATA;
[ x,y ] = fun_Dif_Ur_Modif_Eiler(x0,y0,xn,n);
REPORT;
end % function
function REPORT
global xn n x0 y0 x y;
disp('Diffferencialnoye uravneniye Modif Euler');
disp(' ');
disp('Arguments');
disp(['x0 = ' num2str(x0,'%10.5f') ]);
disp(['y0 = ' num2str(y0,'%10.5f') ]);
disp(['xn = ' num2str(xn,'%10.5f') ]);
disp(['n = ' num2str(n,'%10.5f') ]);
disp('Results');
disp(' i x fx ');
disp(' ______________________________')
i=0;
for i=1:length(x)
xx=x(i);
ffx=y(i);
disp(sprintf('%10.3f\t%10.3f\t %10.3f',i,xx,ffx));
end %for
plot(x,y,'r.');
grid on;
xlabel('x');
ylabel('y');
title('Diffferencialnoye uravneniye Modif Euler');
end % function
Программа. Дифференциальное уравнение.
Метод усовершенствованный Эйлера.
function DATA
global xn n x0 y0;
xn=3;
n=14;
x0=1;
y0=-2;
end
function [ dfdx ] = f(x,y)
dfdx=-y/x;
end
function [ y ] = Usov_Eiler(x,y,h)
arg1=f(x,y);
arg2=f(x+h,y+h*f(x,y));
y=y+h*(arg1+arg2)/2;
end
function [ x,y ] = fun_Dif_Ur_Usov_Eiler(x0,y0,xn,n)
h=(xn-x0)/n;
for i=1:n+1
if i==1
x(i)=x0;
y(i)=y0;
else
x(i)=x(i-1)+h;
y(i)=Usov_Eiler(x(i-1),y(i-1),h);
end %if
end %for i
end % function
function GLAV_Dif_Ur_Usov_Eiler
global xn n x0 y0 x y;
DATA;
[ x,y ] = fun_Dif_Ur_Usov_Eiler(x0,y0,xn,n);
REPORT;
end % function
function REPORT
global xn n x0 y0 x y;
disp('Diffferencialnoye uravneniye Usov Euler');
disp(' ');
disp('Arguments');
disp(['x0 = ' num2str(x0,'%10.5f') ]);
disp(['y0 = ' num2str(y0,'%10.5f') ]);
disp(['xn = ' num2str(xn,'%10.5f') ]);
disp(['n = ' num2str(n,'%10.5f') ]);
disp('Results');
|
|
|
disp(' i x fx ');
disp(' ______________________________')
i=0;
for i=1:length(x)
xx=x(i);
ffx=y(i);
disp(sprintf('%10.3f\t%10.3f\t %10.3f',i,xx,ffx));
end %for
plot(x,y,'r.');
grid on;
xlabel('x');
ylabel('y');
title('Diffferencialnoye uravneniye Usov Euler');
end % function
Программа. Дифференциальное уравнение. Метод Рунге-Кутта.
function DATA
global xn n x0 y0;
xn=3;
n=14;
x0=1;
y0=-2;
end
function [ dfdx ] = f(x,y)
dfdx=-y/x;
end
function [ y ] = Dif_Ur_Runge_Kutt(x,y,h)
k0=f(x,y);
k1=f(x+h/2,y+h*k0/2);
k2=f(x+h/2,y+h*k1/2);
k3=f(x+h,y+h*k2);
y=y+h*(k0+2*k1+2*k2+k3)/6;
end
function [ x,y ] = fun_Dif_Ur_Runge_Kutt(x0,y0,xn,n)
h=(xn-x0)/n;
for i=1:n+1
if i==1
x(i)=x0;
y(i)=y0;
else
x(i)=x(i-1)+h;
y(i)=Dif_Ur_Runge_Kutt(x(i-1),y(i-1),h);
end %if
end %for i
end % function
function GLAV_Dif_Ur_Runge_Kutt
global xn n x0 y0 x y;
DATA;
[ x,y ] = fun_Dif_Ur_Runge_Kutt(x0,y0,xn,n);
REPORT;
end % function
function REPORT
global xn n x0 y0 x y;
disp('Diffferencialnoye uravneniye Runge Kutt');
disp(' ');
disp('Arguments');
disp(['x0 = ' num2str(x0,'%10.5f') ]);
disp(['y0 = ' num2str(y0,'%10.5f') ]);
disp(['xn = ' num2str(xn,'%10.5f') ]);
disp(['n = ' num2str(n,'%10.5f') ]);
disp('Results');
disp(' i x fx ');
disp(' ______________________________')
i=0;
for i=1:length(x)
xx=x(i);
ffx=y(i);
disp(sprintf('%10.3f\t%10.3f\t %10.3f',i,xx,ffx));
end %for
plot(x,y,'r.');
grid on;
xlabel('x');
ylabel('y');
title('Diffferencialnoye uravneniye Runge Kutt');
end % function
Программа. Дифференциальное уравнение.
Сравнение погрешностей методов при одинаковом числе шагов.
function DATA
global xn n x0 y0;
xn=3;
n=4;
x0=1;
y0=-2;
end
function [ dfdx ] = f(x,y)
dfdx=-y/x;
end
function [ y ] = Eiler(x,y,h)
y=y+h*f(x,y);
end
function [ y ] = Modif_Eiler(x,y,h)
xnew=x+h/2;
ynew=y+h*f(x,y)/2;
y=y+h*f(xnew,ynew);
end
function [ y ] = Usov_Eiler(x,y,h)
arg1=f(x,y);
arg2=f(x+h,y+h*f(x,y));
y=y+h*(arg1+arg2)/2;
end
function [ y ] = Dif_Ur_Runge_Kutt(x,y,h)
k0=f(x,y);
k1=f(x+h/2,y+h*k0/2);
k2=f(x+h/2,y+h*k1/2);
k3=f(x+h,y+h*k2);
y=y+h*(k0+2*k1+2*k2+k3)/6;
end
function [ x,yE, yME, yUE, yRK ] = fun_Dif_Ur(x0,y0,xn,n)
h=(xn-x0)/n;
for i=1:n+1
if i==1
x(i)=x0;
yRK(i)=y0;
yE(i)=y0;
yME(i)=y0;
yUE(i)=y0;
else
x(i)=x(i-1)+h;
yRK(i)=Dif_Ur_Runge_Kutt(x(i-1),yRK(i-1),h);
yE(i)=Eiler(x(i-1),yE(i-1),h);
yME(i)=Modif_Eiler(x(i-1),yME(i-1),h);
yUE(i)=Usov_Eiler(x(i-1),yUE(i-1),h);
end %if
end %for i
end % function
function GLAV_Dif_Ur
global xn n x0 y0 x yE yME yUE yRK;
DATA;
[ x,yE, yME, yUE, yRK ] = fun_Dif_Ur(x0,y0,xn,n);
REPORT;
end % function
function REPORT
global xn n x0 y0 x yE yME yUE yRK;
disp('Diffferencialnoye uravneniye ');
disp(' ');
disp('Arguments');
disp(['x0 = ' num2str(x0,'%10.5f') ]);
disp(['y0 = ' num2str(y0,'%10.5f') ]);
disp(['xn = ' num2str(xn,'%10.5f') ]);
disp(['n = ' num2str(n,'%10.5f') ]);
disp('Results');
disp(' i x yE yME yUE yRK ');
disp(' ________________________________________________')
i=0;
for i=1:length(x)
xx=x(i);
ffxE=yE(i);
ffxME=yME(i);
ffxUE=yUE(i);
ffxRK=yRK(i);
disp(sprintf('%10.3f\t%10.3f\t %10.3f\t%10.3f\t %10.3f\t%10.3f',i,xx,ffxE,ffxME,ffxUE,ffxRK));
end %for
plot(x,yE,'r+',x,yME,'b+',x,yUE,'gx',x,yRK,'k+');
grid on;
xlabel('x');
ylabel('y');
title('Diffferencialnoye uravneniye yE r+,yME b+,yUE gx,yRK k+ ');
end % function
Программа. Система дифференциальных уравнений. Метод Рунге-Кутта.
function DATA
global xn n x0 y0;
xn=4;
n=30;
x0=1;
y0=[-2;1;2;-3;3];
end
function [ dfdx ] = f(x,y)
dfdx=[-y(1)/x+y(2); y(2)/x-y(3); y(3)/x; y(4)-y(1)/x; y(5)/x];
end
function [ x,y ] = fun_Dif_Ur_Runge_Kutt(x0,y0,xn,n)
h=(xn-x0)/n
for i=1:n+1
if i==1
x(i)=x0;
y(:,i)=y0;
else
x(i)=x(i-1)+h;
y(:,i)=Dif_Ur_Runge_Kutt(x(i-1),y(:,i-1),h);
end %if
end %for i
end % function
function [ y ] = Dif_Ur_Runge_Kutt(x,y,h)
k0=f(x,y);
k1=f(x+h/2,y+h.*k0/2);
k2=f(x+h/2,y+h.*k1/2);
k3=f(x+h,y+h.*k2);
y=y+h.*(k0+2*k1+2*k2+k3)/6;
end
function GLAV_Dif_Ur_Runge_Kutt
global xn n x0 y0 x y;
DATA;
[ x,y ] = fun_Dif_Ur_Runge_Kutt(x0,y0,xn,n);
REPORT;
end % function
function REPORT
global xn n x0 y0 x y;
disp('Diffferencialnoye uravneniye Runge Kutt');
disp(' ');
disp('Arguments');
disp('______________')
disp('Arguments x0 ');
x0
disp('______________')
disp('Arguments y0 ');
y0
disp('______________')
disp('Arguments xn ');
xn
disp('______________')
disp(['n = ' num2str(n,'%10.5f') ]);
disp('Results');
disp(' x 0 fx ');
disp(' _____________________________________________');
i=0;
lx1=length(x);
z=zeros(lx1,1);
x_fx=[x' z y']
for i=1:length(y(:,1))
switch i
case 1
ssr='-';
case 2
ssr='-.';
case 3
ssr='--';
case 4
ssr='.';
otherwise
ssr='*--';
end%switch
plot(x,y(i,:),ssr);
hold on
end%for i
grid on;
xlabel('x');
ylabel('y ');
title('--1 -.-.2 - - -3...4 - * -5 Diffferencialnoye uravneniye Runge Kutt');
end % function
