Методы решения обыкновенных дифференциальных уравнений

Общий вид обыкновенного дифференциального уравнения, устанавливающего связь между независимой переменной x неизвестной функцией y и ее производными y’,y”,…,y⁽ⁿ⁾, может мыть представлен следующим образом:

Порядок наивысшей производной, входящей в уравнение, называется порядком этого уравнения. Решение дифференциального уравнения (интегрированием) является некоторая функциональная зависимость y=y(x), которая при подстановке в уравнение обращает его в тождество.

Общее решение дифференциального уравнения записывается в виде: y=y(x,c₁,c₂,…,c_n), где c₁,c₂,…,c_n произвольные постоянные.

Решение, полученное из общего решения при фиксированных значениях, называется частным решением уравнения. Постоянные c₁,c₂,…,c_n можно определить, задав n условий. Если эти условия заданы как совокупность значений искомой функции и всех ее производных до (n-1) порядка включительно в некоторой точке x₀, то задача решения уравнения называется задачей Коши, а заданные условия: y(x₀)=y₀, y’(x₀)=y’₀, y”(x₀)=y”₀,…, y^n-1(x₀)=y^n-1₀ называются начальными условия.

Если же условия заданы при нескольких значениях x, то задача решения дифференциального уравнения будет называться граничной или краевой задачей. 

Рассмотрим дифференциальное уравнение первого порядка: соотношение часто удается записать в виде:

Последнее уравнение называется дифференциальным уравнением, разрешенным относительно производной. Значение производной равно тангенсу угла наклона касательной к графику функции в точке (x,y). Функцию f(x,y) будем называть правой частью дифференциального уравнения.

Общим решением уравнения будет являться семейство функций y=y(x,c₁) различающихся значение постоянной c₁. Задаем одно начальное условие y(x₀)=y₀, которое определяет значение c₁и конкретное частное решение – задача Коши.

Для простейшего дифференциального уравнения y’=3x². Общее решение имеет вид y=x³+c, а подставив в общее решение начальное условие x₀=1, y₀=2 вычислим с=1 и определим частное решение как: y=x³+1

Метод Эйлера. Дано дифференциальное уравнение y’=f(x,y), удовлетворяющее начальному условию y(x₀)=y₀. Требуется найти решение на отрезке [a,b]. Разобьем (рис.2.11.1) отрезок интегрирования на n равных частей: x₀=a, x₁= a+h, x₂=x₁+h,…,x_i=x_i–1+h,…,x_n=b, тогда величина шага интегрирования будет равна:

Рис.2.11.1. Метод Эйлера.

Значение функции y₁ в точке x₁ можно определить как точку пересечения касательной проведенной к функции y=y(x) в точке (x₀,y₀) с вертикальной прямой проходящей через точку x₁.  

Тангенс угла наклона касательной есть значение производной в точке (x₀,y₀) и задается правой частью дифференциального уравнения, т.е. tg(β)=f(x₀,y₀). С другой стороны из геометрического представления метода можно записать: Следовательно Откуда тогда и так далее. Решение будет заключаться в последовательном применении формул: где i = 1, 2, 3, …, n Результат будет представлен функцией заданной таблицей.

Модифицированный метод Эйлера. Графическая интерпретация.

Определяем (рис.2.11.2) точку и вычисляем значение функции в этой точке Значение функции y₁ в точке x₁ определяем, как точку пересечения касательной, вычисленной в точке (x_1/2,y_1/2) и проведенной к функции y=y(x) в точке (x₀,y₀), с вертикальной прямой проходящей через точку x₁.

Рис.2.11.2. Модифицированный метод Эйлера.

произвольную точку определим

Система дифференциальных уравнений где x — независимая переменная, а y ₁(x), y ₂(x),..., y_n (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

 

 

 

МАТЛАБ. СБОРНИК ЗАДАЧ