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

Для решения дифференциальных уравнений в форме Коши MatLAB имеет функцию dsolve(‘eqn1’,’eqn2’, …), которая возвращает аналитическое решение системы дифференциальных уравнений с начальными условиями. Они задаются равенствами eqni (вначале задаются уравнения, затем начальные условия).

По умолчанию независимой переменной считается ‘t’ . Можно использовать и другую переменную, включив ее в конец списка параметров функции dsolve. Символ D обозначает производную по независимой переменной, то есть d/dt, при этом D2 означает d^2/dt^2 и т.д.

Начальные условия задаются в виде равенств ‘y(a) = b’ или ‘Dy(a) = b’, где y - независимая переменная, a и b – константы. Если число начальных условий меньше, чем число дифференциальных уравнений, то в решений будут присутствовать произвольные постоянные С1, С2 и т.д. Вывод осуществляется в виде массива записей.

Пример.Решить систему линейных уравнений .

Решение:

>> [x,y]=dsolve('Dx=-x-5*y','Dy=x+y','x(0)=1','y(0)=-1')

>> x = cos(2*t)+2*sin(2*t)

>> y = -cos(2*t)

% построение частного решения

>> t=0:pi/100:2*pi;

>> x=2*sin(2.*t)+cos(2.*t);

>> y=-cos(2.*t);

>> plot(x,y)

>> hold on

>> plot(1,-1,'*')

 

Решение систем дифференциальных уравнений.

 

Для решения систем дифференциальных уравнений существует несколько встроенных процедур. Рассмотрим применение процедуры ode45. Как один из возможных форматов вызова, можно предложить такой: [t,r]=ode45(@DiffEquationFunction,[Tstart,Tfinish], StartVector). Отметим следующее, что процедура ode45 может решить систему уравнений следующего вида: , где - есть векторная функция .

Рассмотрим пример, иллюстрирующий создание исходной функции DiffEquationFunction для вызова ее процедурой ode45. Пусть некоторая точка массы движется в гравитационном поле неподвижных точек с массами и (см. рис.). Уравнение сил в гравитационном поле точек и будет следующим: . Как видим, данное дифференциальное уравнение имеет второй порядок. Но его можно свести к системе дифференциальных уравнений первого порядка: . Будем считать, что данная задача является "плоской" и введем следующие обозначения: , , и . Исходя из таких обозначений, систему дифференциальных уравнений движения точки в гравитационном поле можно представить следующим образом: . Без ущерба для сути решения, значение гравитационной постоянной примем равной 1, , , ,

.

В таком виде систему уравнений можно уже записать как файл-функцию, что мы и сделали, назвав ее threepoint(t,x).

function f=threepoint(t,x)

M1=50; M2=0; C1x=5; C1y=0; C2x=0; C2y=10;

f=[x(3);x(4);...

-M1*(x(1)-C1x)/(sqrt((x(1)-C1x)^2+(x(2)-C1y)^2))^3-...

M2*(x(1)-C2x)/(sqrt((x(1)-C2x)^2+(x(2)-C2y)^2))^3;...

-M1*(x(2)-C1y)/(sqrt((x(1)-C1x)^2+(x(2)-C1y)^2))^3-...

M2*(x(2)-C2y)/(sqrt((x(1)-C2x)^2+(x(2)-C2y)^2))^3];

И решим систему дифференциальных уравнений, вызвав процедуру ode45 из файла-функции dynpoint.m.

function dynpoint()

[t,h]=ode45(@threepoint,[0,1000],[0,0,0,4.3]);

x=h(:,1);

y=h(:,2);

x1=5; y1=0; x2=0; y2=100;

plot(x,y,'b-',x1,y1,'r+',x2,y2,'r*');

При таких начальных параметрах () наша точка движется в поле одного объекта.

dynpoint

То, что мы получили вполне можно назвать приемлемым результатом. Судить о корректности результата с точки зрения физики предоставим физикам. Нашей целью являлось корректная подготовка системы дифференциальных уравнений для ее решения процедурой ode45. С другой стороны, применяя иные встроенные процедуры можно получить отличающися в корне результаты.

Попробуем немного поэкспериментировать и введем значение . Это должно внести возмущения в орбиту движущейся точки.

dynpoint

Крестиком и звездочкой на графике отмечены соответственно и . Помимо процедуры ode45 существует еще ряд встроенных процедур решения систем дифференциальных уравнений. Их описание можно найти в книгах, указанных в списке литературы. Поскольку, решение систем дифференциальных уравнений является очень важным моментом, то мы не ограничимся одним примером.

Рассмотрим следующую задачу, опять же из физики. Рассмотрим траекторию движения пули под действием силы тяжести. При отсутствии сопротивления воздуха, как известно из курса средней школы, это будет парабола. Мы же рассмотрим случай, когда сила сопротивления воздуха пропорциональна квадрату скорости и противоположна направлению движения. Как говорит школьный курс физики, уравнение баланса сил будет следующим: . Учитывая то, что ускорение – производная скорости по времени, распишем это уравнение в векторном виде: . Где - плотность воздуха, - масса пули, - площадь поперечного сечения. Распишем это уравнение покоординатно: . В таком виде система дифференциальных уравнений готова для того, чтобы попытаться решить ее при помощи процедуры ode45. Пусть масса пули – 10 грамм, поперечник – 1 см, плотность воздуха – 1 килограмм на кубический метр. С такими данными мы составим файл-функцию airpoint.m.

function u=airpoint(t,v)

g=10; ro=1; s=0.0001; m=0.01; k=ro*s/2/m;

u=[-k*sqrt(v(1)^2+v(2)^2)*v(1);-g-k*sqrt(v(1)^2+v(2)^2)*v(2)];

В такой системе уравнений аргументом являются скорость и время. Если мы хотим найти траекторию движения, то мы должны принять во внимание: . Полученные массивы точек , , можно в дальнейшем обработать. Произвести интерполяцию, аппроксимацию и т.д. Но мы интегрирование заменим суммированием: . В этом случае начальный момент времени равен 0, пуля находится в начале координат. Создадим файл-функцию dynairpoint, в которой вызовем процедуру ode45. В начальный момент времени скорость пули по горизонтали 800, по вертикали – 100.

function dynairpoint()

[t,h]=ode113(@airpoint,[0,5.6],[800,100]);

vx=h(:,1); vy=h(:,2);

m=length(t);

x0=0; y0=0;

x(1)=x0; y(1)=0;

for i=2:m

x(i)=x(i-1)+vx(i-1)*(t(i)-t(i-1));

y(i)=y(i-1)+vy(i-1)*(t(i)-t(i-1));

end;

[t1,h1]=ode113(@airpoint,[0,5.6],[800,100]);

vx1=h1(:,1); vy1=h1(:,2);

m1=length(t1);

x01=0; y01=0;

x1(1)=x0; y1(1)=0;

for i=2:m1

x1(i)=x1(i-1)+vx1(i-1)*(t1(i)-t1(i-1));

y1(i)=y1(i-1)+vy1(i-1)*(t1(i)-t1(i-1));

end;

figure

plot(x,y,'r+',x1,y1)

grid on

Следует отметить, что диапазон времени [0;5.6] был эмпирическим путем подобран так, что траектория движения отображена от момента вылета до падения. Уровень "земли" соответствует значению ординаты 0. Момент прикосновения к "земле" можно и подсчитать. Оставим это на Ваше усмотрение. В m-файле dynairpoint.m вызвана не только процедура ode45, но и ode113. И Вы можете сравнить результаты.

dynairpoint

 

 

Мы решали систему дифференциальных уравнений . Но можно решить и такую систему уравнений: . Произведем небольшие изменения и напишем файл-функцию airpoint1.m.

function u=airpoint1(t,x)

g=10; ro=1; s=0.0001; m=0.01; k=ro*s/2/m;

u=[x(3);x(4);...

-k*sqrt(x(3)^2+x(4)^2)*x(3);...

-g-k*sqrt(x(3)^2+x(4)^2)*x(4)];

Для решения системы четырех дифференциальных уравнений, которая записана файл-функцией airpoint1.m создадим файл-функцию dynairpoint1.m.

function dynairpoin1t()

[t,h]=ode45(@airpoint1,[0,5.6],[0,0,800,100]);

x=h(:,1); y=h(:,2);

[t1,h1]=ode113(@airpoint1,[0,5.6],[0,0,800,100]);

x1=h1(:,1); y1=h1(:,2);

figure

plot(x,y,'r-',x1,y1,'b-')

grid on

dynairpoint1

В этом случае практически нет различий между решениями ode45 и ode113.