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

   1. Пример 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.0)

Для этого решаем дифференциальное уравнение Эйлера относительно 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.0)

где в нашем случае 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' ) % заголовок