Примеры выполнения заданий

A. Пример 11a

Решить методом начальных параметров пример 1a. Сравнить решение с аналитическим. Построить графики.

Используем программу для примера 1a для составления этой программы. Вводим исходные данные. Решаем пример 1a и заполняем таблицу.

clear all

format long

disp('Решаем пример 11a')

nnp = 10; % число точек интегрирования

syms x y Dy D2y % описали символические переменные

F = x^2+y^2+Dy^2; % подвнтегральная функция

x1 = -1;

y1 = 1;

x2 = 1;

y2 = 2;

fprintf('Подынтегральная функция: F=%s\n',char(F))

fprintf('Граничные условия: y(%d)=%d; y(%d)=%d\n',x1,y1,x2,y2)

dFdy = diff(F,y);

dFdy1 = diff(F,Dy);

d_dFdy1_dx = diff(dFdy1,x);

d_dFdy1_dy = diff(dFdy1,y);

d_dFdy1_dy1 = diff(dFdy1,Dy); % d(dF/dy')/dy'

dFy1dx = d_dFdy1_dx + d_dFdy1_dy*Dy + d_dFdy1_dy1*D2y;

Euler = simple(dFdy-dFy1dx);

deqEuler = [ char(Euler) '=0' ]; % составили уравнение

Sol = dsolve(deqEuler,'x'); % решаем уравнение Эйлера

if length(Sol)~=1 % решений нет или более одного

error('Нет решений или более одного решения!');

end

SolLeft = subs(Sol,x,sym(x1));

SolRight = subs(Sol,x,sym(x2));

EqLeft = [char(SolLeft) '=' char(sym(y1))];

EqRight = [char(SolRight) '=' char(sym(y2))];

Con = solve(EqLeft,EqRight); % решаем систему

C1 = Con.C1;

C2 = Con.C2;

Sol1a = vpa(eval(Sol),14);

xpl = linspace(x1,x2);

y1a=subs(Sol1a,x,xpl);

Решаем пример 11a

Подынтегральная функция: F=x^2+y^2+Dy^2

Граничные условия: y(-1)=1; y(1)=2

Для применения решателей систем дифференциальных уравнений приводим уравнение Эйлера 2-го порядка к системе 2-х дифференциальных уравнений 1-го порядка путём замены

(11.4)

Для этого решаем дифференциальное уравнение Эйлера относительно y¢¢ и формируем правые части системы дифференциальных уравнений. Записываем их в файл.

f2 = solve(deqEuler,D2y); % решаем относительно y''

f2 = subs ( f2, {y,Dy}, {sym('y(1)'),sym('y(2)')} );

rp{1} = 'function dydx = MyRightPart(x,y)';

rp{2} = 'dydx=zeros(2,1);';

rp{3} = 'dydx(1)=y(2);';

rp{4} = [ 'dydx(2)=' char(f2) ';' ];

disp('Текст файла MyRightPart.m')

fprintf ( '%s\n', rp{:} );

fid = fopen ( 'C:\Iglin\Matlab\MyRightPart.m', 'w' );

fprintf ( fid, '%s\n', rp{:} );

fclose(fid); % закрываем файл

Текст файла MyRightPart.m

function dydx = MyRightPart(x,y)

dydx=zeros(2,1);

dydx(1)=y(2);

dydx(2)=y(1);

Мы сформировали систему дифференциальных уравнений

(11.5)

где в нашем случае x₁=-1; x₂=1. Если бы было известно начальное условие y₂(x₁), то эту систему можно было бы решить с помощью стандартных численных методов. Обозначим y₂(x₁) как неизвестную величину: t=y₂(x₁). Присвоив ей какое-либо пробное значение, можно решить систему дифференциальных уравнений и найти функцию f=y₁(x₂)-2. Очевидно, f можно рассматривать как функцию от t. То есть нужно решить уравнение f(t)=0. В нашем случае, когда исходная система дифференциальных уравнений является линейной, уравнение относительно t также будет линейным. То есть функция f(t) имеет структуру f(t)=at+b. Чтобы построить эту функцию, нужно решить 2 начальные задачи для t=0 и t=1. Решаем эти задачи.

xr = linspace(x1,x2,nnp+1); % точки для численного решения

y0 = [y1,0]; % решаем СДУ при y'(x1)=0;

[xx,YY] = ode45('MyRightPart',xr,y0);

yend0 = YY(nnp+1,1)-y2

y0 = [y1,1]; % решаем СДУ при y'(x1)=1;

[xx,YY] = ode45('MyRightPart',xr,y0);

yend1 = YY(nnp+1,1)-y2

yend0 =

1.76219612254200

yend1 =

5.38905701685360

Система линейных уравнений (11.3) в данном случае состоит из одного уравнения. Для нахождения неизвестного t=y₂(x₁) проводим линейную интерполяцию.

y0 = [y1,yend0/(yend0-yend1)]

y0 =

1.00000000000000 -0.48587364497652

Решаем систему дифференциальных уравнений при найденных действительных начальных условиях. Строим график.

[xx,YY] = ode45('MyRightPart',xr,y0);

plot ( xpl,y1a,'--b', xr,YY(:,1),'-r' )

title ( '\bfExample 11a' ) % заголовок

xlabel( 'x')

ylabel('y(x)') % метки осей

Рис. 11.2. Решение примера 11a

Ответ. График экстремали показан на рис.11.2 сплошной красной линией. Штриховая синяя линия – решение примера 1a. Видно, что в точках, где печатается численное решение, оно сливается с аналитическим. Неизвестное начальное условие: y¢(x₁)=-0.48587364.

B. Пример 11b

Решить методом начальных параметров пример 2. Сравнить решение с аналитическим. Построить графики.

Программу для этого примера напишем на основе программы для примера 2 с использованием программы для примера 11a. Находим аналитическое решение примера 2. Заполняем таблицу.

clear all

format long

disp('Решаем пример 11b')

nnp = 10; % число точек интегрирования

syms x y z Dy D2y Dz D2z % описали переменные

F = Dy^2+Dz^2+2*y*z; % подынтегральная функция

x1 = -2;

y1 = 1;

z1 = 0;

x2 = 2;

y2 = 0;

z2 = 2;

fprintf('Подынтегральная функция: F=%s\n',char(F))

fprintf('Граничные условия слева: y(%d)=%d; z(%d)=%d\n',x1,y1,x1,z1)

fprintf('Граничные условия справа: y(%d)=%d; z(%d)=%d\n',x2,y2,x2,z2)

dFdy = diff(F,y);

dFdy1 = diff(F,Dy);

d_dFdy1_dx = diff(dFdy1,x); % d(dF/dy')/dx

d_dFdy1_dy = diff(dFdy1,y);

d_dFdy1_dy1 = diff(dFdy1,Dy);

d_dFdy1_dz = diff(dFdy1,z);

d_dFdy1_dz1 = diff(dFdy1,Dz);

dFy1dx = d_dFdy1_dx + d_dFdy1_dy*Dy + d_dFdy1_dy1*D2y + d_dFdy1_dz*Dz + d_dFdy1_dz1*D2z;

dFdz = diff(F,z);

dFdz1 = diff(F,Dz);

d_dFdz1_dx = diff(dFdz1,x); % d(dF/dz')/dx

d_dFdz1_dy = diff(dFdz1,y);

d_dFdz1_dy1 = diff(dFdz1,Dy);

d_dFdz1_dz = diff(dFdz1,z);

d_dFdz1_dz1 = diff(dFdz1,Dz);

dFz1dx = d_dFdz1_dx + d_dFdz1_dy*Dy + d_dFdz1_dy1*D2y + d_dFdz1_dz*Dz + d_dFdz1_dz1*D2z;

EulerY = simple(dFdy-dFy1dx);

EulerZ = simple(dFdz-dFz1dx);

deqEulerY = [char(EulerY) '=0']; % уравнение Y

deqEulerZ = [char(EulerZ) '=0']; % уравнение Z

Sol = dsolve(deqEulerY,deqEulerZ,'x'); % решаем

if length(Sol)~=1 % решений нет или более одного

error('Нет решений или более одного решения!');

end

SolLeftY = subs(Sol.y,x,sym(x1)); % x1 в y

SolLeftZ = subs(Sol.z,x,sym(x1)); % x1 в z

SolRightY = subs(Sol.y,x,sym(x2)); % x2 в y

SolRightZ = subs(Sol.z,x,sym(x2)); % x2 в z

EqLeftY = [char(vpa(SolLeftY,14)) '=' char(sym(y1))];

EqLeftZ = [char(vpa(SolLeftZ,14)) '=' char(sym(z1))];

EqRightY = [char(vpa(SolRightY,14)) '=' char(sym(y2))];

EqRightZ = [char(vpa(SolRightZ,14)) '=' char(sym(z2))];

Con = solve(EqLeftY,EqLeftZ,EqRightY,EqRightZ);

C1 = Con.C1;

C2 = Con.C2;

C3 = Con.C3;

C4 = Con.C4;

Sol2Y = vpa(eval(Sol.y),14);

Sol2Z = vpa(eval(Sol.z),14);

xpl = linspace(x1,x2); % массив абсцисс

y2a = subs(Sol2Y,x,xpl);

z2a = subs(Sol2Z,x,xpl);

Решаем пример 11b

Подынтегральная функция: F=Dy^2+Dz^2+2*y*z

Граничные условия слева: y(-2)=1; z(-2)=0

Граничные условия справа: y(2)=0; z(2)=2

Сведём систему 2-х дифференциальных уравнений Эйлера 2-го порядка к системе 4-х нормальных дифференциальных уравнений 1-го порядка вида (11.1). Решим уравнения Эйлера относительно y¢¢, z¢¢. Подставим в оба уравнения y, z, y¢, z¢. Сформируем правые части для системы дифференциальных уравнений и запишем их в файл.

f2yz = solve(deqEulerY,deqEulerZ,D2y,D2z);

f2=subs(f2yz.D2y,{y,Dy,z,Dz},{sym('y(1)'),sym('y(2)'),sym('y(3)'),sym('y(4)')});

f4=subs(f2yz.D2z,{y,Dy,z,Dz},{sym('y(1)'),sym('y(2)'),sym('y(3)'),sym('y(4)')});

rp{1} = 'function dydx = MyRightPart(x,y)';

rp{2} = 'dydx=zeros(4,1);';

rp{3} = 'dydx(1)=y(2);';

rp{4} = [ 'dydx(2)=' char(f2) ';' ];

rp{5} = 'dydx(3)=y(4);';

rp{6} = [ 'dydx(4)=' char(f4) ';' ];

disp('Текст файла MyRightPart.m')

fprintf('%s\n',rp{:});

fid = fopen ( 'C:\Iglin\Matlab\MyRightPart.m', 'w' );

fprintf(fid,'%s\n',rp{:});

fclose(fid); % закрываем файл

Текст файла MyRightPart.m

function dydx = MyRightPart(x,y)

dydx=zeros(4,1);

dydx(1)=y(2);

dydx(2)=y(3);

dydx(3)=y(4);

dydx(4)=y(1);

Применим метод начальных параметров для решения задачи. Неизвестные у нас обозначены

(11.6)

В начальной точке x₁ неизвестны y₂(x₁)=t₁ и y₄(x₁)=t₂. Найдём их из решения системы 2-х уравнений

(11.7)

Система дифференциальных уравнений Эйлера у нас является линейной, поэтому и система уравнений (11.7) также будет линейной.

(11.8)

Найдём коэффициенты этой системы и правые части.

xr = linspace(x1,x2,nnp+1); % точки для численного решения

A = zeros(2,2); % матрица системы уравнений

b = zeros(2,1); % вектор правых частей

y0 = [y1;0;z1;0]; % решаем СДУ при y'(x1)=0; z'(x1)=0;

[xx,YY] = ode45('MyRightPart',xr,y0);

b=YY(nnp+1,[1 3])' - [y2;z2]

y0 = [y1;1;z1;0]; % решаем СДУ при y'(x1)=1; z'(x1)=0;

[xx,YY] = ode45('MyRightPart',xr,y0);

A(:,1) = YY(nnp+1,[1 3])'-[y2;z2]-b;

y0 = [y1;0;z1;1]; % решаем СДУ при y'(x1)=0; z'(x1)=1;

[xx,YY] = ode45('MyRightPart',xr,y0);

A(:,2) = YY(nnp+1,[1 3])'-[y2;z2]-b

b =

13.32736606398895

11.98099902570117

A =

13.26662430099436 14.02342479230629

14.02342479230629 13.26662430099436

Решаем систему линейных уравнений (11.8), находим недостающие начальные условия. Для этих начальных условий решаем систему дифференциальных уравнений вида (11.1) и строим график полученного решения.

yz0 = -A\b; % нашли начальные условия

y0 = [y1;yz0(1);z1;yz0(2)] % истинные начальные условия

[xx,YY] = ode45('MyRightPart',xr,y0); % решаем

plot3(xpl,y2a,z2a,'--b',xr,YY(:,1),YY(:,3),'-r')

title ( '\bfExample 11b' ) % заголовок

xlabel('x')

ylabel('y(x)')

zlabel('z(x)')

view(205,30)

grid on

box on

y0 =

1.00000000000000

0.42582034231673

-1.35320471612876

Рис. 11.3. Решение примера 11b

Ответ. График экстремали показан на рис.11.3 сплошной красной линией. Он практически сливается с решением примера 2, которое показано штриховой синей линией. Неизвестные начальные условия: y¢(x₁)=0.42582; z¢(x₁)=-1.35320.

C. Пример 11c

Решить методом начальных параметров пример 3. Сравнить решение с аналитическим. Построить графики.

Уравнение Эйлера (3.6) является уравнением 4-го порядка с 2-мя граничными условиями на левом конце и 2-мя на правом. Сведём его к нормальной системе 4-х уравнений 1-го порядка заменой

(11.9)

Неизвестные начальные условия t₁=y₃(x₁) и t₂=y₄(x₁) найдём из решения системы уравнений

(11.10)

Так как уравнение (3.6) линейное, то и система (11.10) будет линейной. Решая её, найдём начальные условия, а затем и решение уравнения Эйлера.

Для составления программы используем программы для примеров 3 и 11b. Вначале повторяем решение примера 3.

clear all

format long

disp('Решаем пример 11c')

nnp = 10;

syms x y Dy D2y D3y D4y % описали переменные

F = D2y^2-2*Dy^2+4*y*Dy+y^2-2*y*sin(x);

x1 = -1;

y1 = 1;

Dy1 = -1;

x2 = 1;

y2 = 2;

Dy2 = 1;

fprintf('Подынтегральная функция: F=%s\n',char(F))

fprintf('Граничные условия слева: y(%d)=%d; y''(%d)=%d\n',x1,y1,x1,Dy1)

fprintf('Граничные условия справа: y(%d)=%d; y''(%d)=%d\n',x2,y2,x2,Dy2)

dFdy = diff(F,y);

dFdy1 = diff(F,Dy);

dFdy2 = diff(F,D2y); % dF/dy''

d_dFdy1_dx = diff(dFdy1,x); % d(dF/dy')/dx

d_dFdy1_dy = diff(dFdy1,y); % d(dF/dy')/dy

d_dFdy1_dy1 = diff(dFdy1,Dy); % d(dF/dy')/dy'

d_dFdy1_dy2 = diff(dFdy1,D2y); % d(dF/dy')/dy''

dFy1dx = d_dFdy1_dx + d_dFdy1_dy * Dy + d_dFdy1_dy1 * D2y + d_dFdy1_dy2*D3y;

d_dFdy2_dx = diff(dFdy2,x); % d(dF/dy'')/dx

d_dFdy2_dy = diff(dFdy2,y); % d(dF/dy'')/dy

d_dFdy2_dy1 = diff(dFdy2,Dy); % d(dF/dy'')/dy'

d_dFdy2_dy2 = diff(dFdy2,D2y); % d(dF/dy'')/dy''

dFy2dx = d_dFdy2_dx + d_dFdy2_dy * Dy + d_dFdy2_dy1 * D2y + d_dFdy2_dy2 * D3y;

d_dFdy2dx_dx = diff(dFy2dx,x); % d((dFy'')/dx)/dx

d_dFdy2dx_dy = diff(dFy2dx,y); % d((dFy'')/dx)/dy

d_dFdy2dx_dy1 = diff(dFy2dx,Dy); % d((dFy'')/dx)/dy'

d_dFdy2dx_dy2 = diff(dFy2dx,D2y); % d((dFy'')/dx)/dy''

d_dFdy2dx_dy3 = diff(dFy2dx,D3y); % d((dFy'')/dx)/dy'''

d2Fy2dx2 = d_dFdy2dx_dx + d_dFdy2dx_dy * Dy + d_dFdy2dx_dy1 * D2y + d_dFdy2dx_dy2 * D3y + d_dFdy2dx_dy3 * D4y;

Euler = simple(dFdy-dFy1dx+d2Fy2dx2);

deqEuler = [char(Euler) '=0']; % составили уравнение

Sol = dsolve ( deqEuler, 'x' );

if length(Sol)~=1 % решений нет или более одного

error('Нет решений или более одного решения!');

end

dydx = diff(Sol,x); % нашли производную

slY = subs(Sol,x,sym(x1));

slDY = subs(dydx,x,sym(x1));

srY = subs(Sol,x,sym(x2));

srDY = subs(dydx,x,sym(x2));

elY = [char(vpa(slY,14)) '=' char(sym(y1))];

elDY = [char(vpa(slDY,14)) '=' char(sym(Dy1))];

erY = [char(vpa(srY,14)) '=' char(sym(y2))];

erDY = [char(vpa(srDY,14)) '=' char(sym(Dy2))];

Con = solve(elY,elDY,erY,erDY);

C1 = Con.C1;

C2 = Con.C2;

C3 = Con.C3;

C4=Con.C4;

Sol3 = vpa(eval(Sol),14); % подставляем C1-C4;

xpl = linspace(x1,x2);

y3 = subs(Sol3,x,xpl);

Решаем пример 11c

Подынтегральная функция: F=D2y^2-2*Dy^2+4*y*Dy+y^2-2*y*sin(x)

Граничные условия слева: y(-1)=1; y'(-1)=-1

Граничные условия справа: y(1)=2; y'(1)=1

Разрешаем уравнение Эйлера относительно y^IV и подставляем в него обозначения (11.9). Формируем правые части для системы дифференциальных уравнений, к которой сводится уравнение Эйлера. Записываем их в файл.

f4 = solve(deqEuler,D4y); % находим D4y

f4 = subs(f4,{y,Dy,D2y,D3y},{sym('y(1)'),sym('y(2)'),sym('y(3)'),sym('y(4)')});

rp{1} = 'function dydx = MyRightPart(x,y)';

rp{2} = 'dydx=zeros(4,1);';

rp{3} = 'dydx(1)=y(2);';

rp{4} = 'dydx(2)=y(3);';

rp{5} = 'dydx(3)=y(4);';

rp{6} = [ 'dydx(4)=' char(f4) ';' ];

disp('Текст файла MyRightPart.m')

fprintf('%s\n',rp{:});

fid = fopen ( 'C:\Iglin\Matlab\MyRightPart.m', 'w' );

fprintf(fid,'%s\n',rp{:});

fclose(fid); % закрываем файл

Текст файла MyRightPart.m

function dydx = MyRightPart(x,y)

dydx=zeros(4,1);

dydx(1)=y(2);

dydx(2)=y(3);

dydx(3)=y(4);

dydx(4)=-y(1)+sin(x)-2*y(3);

Формируем коэффициенты и свободные члены системы линейных уравнений (11.3, 11.10). Для этого 3 раза решаем начальную задачу при значениях неизвестных начальных параметров: {t₁,t₂}={0,0}; {t₁,t₂}={1,0}; и {t₁,t₂}={0,1}. Решаем систему (11.10) – находим неизвестные начальные параметры. Решаем начальную задачу при этих значениях начальных параметров. Рисуем график.

xr = linspace(x1,x2,nnp+1); % точки для численного решения

A = zeros(2,2);

b = zeros(2,1);

y0 = [y1;Dy1;0;0]; % начальные условия (0,0)

[xx,YY] = ode45('MyRightPart',xr,y0); % решаем СДУ

b = YY(nnp+1,[1 2])' - [y2;Dy2] % правые части

y0 = [y1;Dy1;1;0]; % начальные условия (1,0)

[xx,YY] = ode45('MyRightPart',xr,y0); % решаем СДУ

A(:,1) = YY(nnp+1,[1 2])'-[y2;Dy2]-b; % 1-й столбец матрицы A

y0 = [y1;Dy1;0;1]; % начальные условия (0,1)

[xx,YY] = ode45('MyRightPart',xr,y0); % решаем СДУ

A(:,2) = YY(nnp+1,[1 2])'-[y2;Dy2]-b % 2-й столбец матрицы A

yz0 = -A\b; % нашли начальные условия

y0 = [y1;Dy1;yz0] % истинные начальные условия

[xx,YY] = ode45('MyRightPart',xr,y0); % решаем

plot(xpl,y3,'--b',xr,YY(:,1),'-r' )

title ( '\bfExample 11c' ) % заголовок

xlabel('x')

ylabel('y(x)')

b =

-3.54609584741614

-2.69446117363805

A =

0.90929722937057 0.87079594142352

0.03850131148804 0.90929725087375

y0 =

1.00000000000000

-1.00000000000000

1.10693966258040

2.91636485466376

Рис. 11.4. Решение примера 11c

Ответ. График экстремали показан на рис.11.4 сплошной красной линией. Он практически сливается с решением примера 3, которое показано штриховой синей линией. Неизвестные начальные условия: y¢¢(x₁)=1.10694; y¢¢¢(x₁)=2.91636.

Задание

Решить примеры 1a , 2 и 3 методом начальных параметров. Сравнить решения с аналитическими.