Обыкновенные дифференциальные уравнения и системы.

Дифференциальным уравнением n-го порядка называется уравнение: H(t, x, x^’, x^’’, …, x⁽ⁿ⁾)=0

Решением дифференциального уравнения называется функция x(t), которая обращает уравнение в тождество.

Системой дифференциальных уравнений n-го порядка называется система вида:

x1^’=f1(t, x1, x2, x3, …, xn) x2^’=f1(t, x1, x2, x3, …, xn) ………………………….. xn^’=f1(t, x1, x2, x3, …, xn)

Системой линейных дифференциальных уравнений называют систему вида:

x1^’= xj + b1 x2^’= xj + b1 .………………….. xn^’= xj + b1

Решением системы является вектор x(t)=x₁(t), x₂(t), …,x_n(t), который обращает вышеприведенные системы уравнений в тождества.

Различают: задачу Коши, когда все дополнительные условия заданы в одной (чаще всего в начальной) точке интервала значений t; краевую задачу, когда дополнительные условия заданы на границах интервала значений t.

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

· f вектор-функция для вычисления правой части уравнения вида x^’=f(x,t), удовлетворяющей начальному условию x(t₀)=x₀ или системы дифференциальных уравнений, данных в матричном виде dx/dt=f(t,x) с начальным условием x(t₀)=x₀ данным в матричном виде.

· interval массив из двух чисел, определяющий интервал интегрирования дифференциального уравнения или системы.

· X0 вектор начальных условий системы дифференциальных уравнений.

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

Все функции возвращают:

· Массив Т содержащий координаты узлов сетки, в которых ищется решение.

· Матрицу Х в которой i-й столбец является значением вектор-функции решения в узле T_i.

 

Решим задачу Коши: dx/dt + x =sin xt при  x(0)=1.5. Напишем М-функцию. Используем функцию ode113, которая реализует метод Адамса. Она строит график (рис.4.1) решения.

Рис.4.1. График решения дифференциального уравнения

function y=D(t,x) y=-x+sin(x*t); end

В качестве ее параметров используются: @D ссылка на М-функцию, [0 35] интервал интегрирования уравнения, начальное условие x(0)=1.5 указано в виде 1,5 на соответствующем месте в списке параметров. Обращение к функции: ode113(@D, [0 35], 1.5). 

Решим систему дифференциальных уравнений на интервале [0 10] при

y1^/=cos(y1y2)       y2^/=sin(y1+ty2)

y1(0)=0 y2(0)=0

Напишем М-функцию syst.

function dy=syst(t,y) dy=zeros(2,1); % zeros это матрица из нулей, указанной в % скобках размерности dy(1)=cos(y(1)*y(2)); dy(2)=sin(y(1)+y(2)*t); end

Используем функцию ode23 решающую методом Рунге-Кутта 2-3 порядка. Ее вызов [T, Y]=ode23(@syst, [0 10], [0 0]); где явно указаны выходные параметры вектор Т и матрица Y, а также указаны интервал [0 10] и начальные значения [0 0].

Затем формируем график решения plot(T,Y(:,1)’-r’,T,Y(:,2),’-k‘). Получаем (рис.4.2) график.

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

Решение жестких задач. Для решения жестких задач действуют в два этапа. Сначала убеждаются, что система жесткая (действительные части всех собственных чисел матрицы dx / dt = B отрицательны и величина s = max | Re (h)|/ min | Re (h)|, называемая числом жесткости системы, велика). После этого решают систему с помощью соответствующих функций матпакета.

 

Найдем решение задачи Коши для следующей системы:

 

119.46   185.38     126.88       121.03 dX = -10.395  -10.136 -3.636    8.577    X dt       -53.302 -85.932 -63.182  54.211                 -115.58 -181.75  -112.8    -199

имеющей начальное условие X(0)=(1 1 1 1);

Опишем систему в функции exSyst_1.

function dx=exSyst_1(t,x) B=[119.46, 185.38, 126.88, 121.03; -10.395, -10.136, -3.636, 8.577; -53.302, -85.932, -63.182, 54.211; -115.58, -181.75, -112.8, -199]; dx=B*x; end

Сначала проверим, является ли матрица В матрицей жесткой системы. Оператор eig(B) возвращает собственные числа матрицы. Одно из них неотрицательное, то есть условие соблюдено не полностью.

Далее находим число жесткости s=28, что не очень много, то есть система не очень жесткая.

Затем команда ode15s решает задачу. Возвращается вектор Т и матрица Х, число столбцов в которой равно числу переменных.

Строим график (рис.4.3). На нем видно, что корни изменяются в зависимости от параметра Т как бы синусоидально. Такой характер изменений корней характерен для жестких систем. По мере роста числа жесткости, корни колеблются (в зависимости от параметра Т) все более резко, образуя «детский рисунок» на экране.

function Start_difur B=[119.46, 185.38, 126.88, 121.03; -10.395, -10.136, -3.636, 8.577; -53.302, -85.932, -63.182, 54.211; -115.58, -181.75, -112.8, -199] L=eig(B) s=max(real(abs(L)))/min(real(abs(L))) [T, X]=ode15s(@exSyst_1,[0 0.1],[1;1;1;1]); plot(T,X(:,1),'r-',T,X(:,2),'g-',T,X(:,3),'c-',T,X(:,4),'k-'); end B = 119.4600 185.3800 126.8800 121.0300 -10.3950 -10.1360 -3.6360 8.5770 -53.3020 -85.9320 -63.1820 54.2110  -115.5800 -181.7500 -112.8000 -199.0000 >> L=eig(B) L =  -64.0747 +74.6101i  -64.0747 -74.6101i  -28.2167          3.5081          s = 28.0342

Рис.4.3. График решения жесткой системы.

Решение дифференциальных уравнений высших порядков. Если в уравнении есть производные высших порядков, то действуем, преобразуя уравнение с понижением степени производной, но с увеличением числа переменных. Решим краевую задачу. d ² x / dt ² +4 dx / dt +13 x = e^{sin t}

при x (0.25)=-1, x ^/ (0.25)=1 на интервале [0.25; 2].

Преобразуем следующим образом. Пусть y=dx/dt. Тогда получим систему уравнений. 

dy/dt=-4y-13x+e^sin t dx/dt=y

с начальными условиями y(0.25)=1, x(0.25)=-1;

опишем ее в функции exSistVP.

function f=exSistVP(t,x) f=[-4*x(1)-13*x(2)+exp(sin(t)); x(1)]; end

Для ее решения воспользуемся функцией ode45.

function Start_difur2 x0=[1;-1];  [T, X]=ode45(@exSistVP,[0.25 2], x0); plot(T, X(:,1),'r-',T,X(:,2),'g-') end

Ей передадим имя М-функции для вышеописанной системы, интервал изменения параметра Т, вектор х0 начальных условий для переменных

Функция ode45 вернет нам матрицу для Х и вектор для Т. Построим (рис.4.4) график.

 

Рис.4.4. Решение диф. уравнения высших порядков.