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

Дифференциальным уравнением 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, которая реализует метод Адамса. Она строит график (рис.1.30) решения.

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

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‘). Получаем (рис.1.31) график.

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

Решение жестких задач. Для решения жестких задач действуют в два этапа. Сначала убеждаются, что система жесткая (действительные части всех собственных чисел матрицы 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 решает задачу. Возвращается вектор Т и матрица Х, число столбцов в которой равно числу переменных.

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

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

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

Решение дифференциальных уравнений высших порядков. Если в уравнении есть производные высших порядков, то действуем, преобразуя уравнение с понижением степени производной, но с увеличением числа переменных. Решим краевую задачу. 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 вернет нам матрицу для Х и вектор для Т. Построим (рис.1.33) график.

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

Обработка результатов эксперимента. Построение полиномиальной зависимости. Для этого есть функция polyfit, которой передаются массив абсцисс, массив ординат (экспериментальных точек) и степень полинома, который требуется получить так, чтобы он описывал экспериментальную зависимость массива ординат от массива абсцисс. Функция возвращает вектор коэффициентов вышеупомянутого полинома.

function Start_poly10 x=[0 4 10 15 21 29 36 51 68]; y=[66.7 71.0 76.3 80.6 85.7 92.9 99.4 113.6 125.1]; [a]=polyfit(x,y,1) xmas=0:0.1:68; ymas=a(1)*xmas+a(2); % для графика plot(x,y,'o',xmas,ymas,'-r'); end

%коэффициенты % полинома a = 0.8706 67.5078

Будет построен график (рис.1.34),на котором можно будет сравнить исходные точки и расчетный график. Расчетный график построен по коэффициентам полинома.

Рис.1.34. Данные эксперимента и расчетный график.

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

Если же требуется подобрать такую зависимость, которая не содержит некоторые степени (например, y =3 x ⁵ +6 x ³ + x +65) то тем самым требуется найти минимум функции. Поясним это на примере. Пусть требуется функция вида y = c ₁ + c ₂ x ² + c ₃ x ³, тогда, используя функцию fminsearch (‘ fname ’, c) получим вектор коэффициентов.

Однако в функцию fminsearch следует передавать не вышеописанную функцию у(х), а переделанный ее вариант вида y_i^experimental -(c ₁ + c ₂ x_i ² + c ₃ x_i ³). Тогда минимум функции будет достигнут в случае, если вычисленные ее значения практически совпадают с экспериментальными, то есть идеально описывают экспериментальные данные. Опишем для этого функцию exGl.

function s=exGl(c) global x; global y; s=0; for i=1:length(x)    s=s+(y(i)-c(1)-c(2)*x(i)^2-c(3)*x(i)^3)^2; end end
global x; global y; c=[2, 1, 0];	Сначала определяем вектор предполагаемых приблизительно значений с (если его неудачно выбрать то ответ не получится). Потом объявляем переменные х и у глобальными.

Затем вводим экспериментальные данные. После этого запускаем функцию fminsearch. Полученный от нее вектор значений с используем для построения (рис.1.35) расчетной кривой. Потом сравниваем, насколько точно кривая идет через экспериментальные точки.

Рис.1.34. Расчетная кривая идет между экспериментальными точками.

function Start_3 x=[10.1 10.2 10.3 10.8 10.9 11 11.1 11.4 12.2 13.3 13.8 14 14.4 14.5 15 15.6 15.8 17 18.1 19]; y=[24 36 26 45 34 37 55 51 75 84 74 91 85 87 94 92 96 97 98 99];
[c,err]=fminsearch('exGl',c) t=x(1):0.1:x(length(x)); fy=c(3)*t.^3+c(2)*t.^2+c(1); plot(x,y,'o',t,fy,'r-'); end	c = -93.7295 1.9899 -0.0774 err = 754.2148

Интерполяция функций. Для реализации сплайн-интерполяции используется функция yi = interp 1(x, y, xi, method) где х массив абсцисс, у массив ординат (экспериментальных точек), xi точки, в которых необходимо вычислить значение с помощью сплайна, yi функция в этих точках. В зависимости от используемых методов, получаются различные результаты (Рис.1.36). Для уплотнения экспериментальных точек, создан массив xi из точек с одинаковым шагом, и для них будут найдены значения функции y(x). Кроме того, создан массив из точек t и эти точки расположены между экспериментальными. Для них тоже будут найдены значения функции y(t)

Рис.1.36. В зависимости от используемых методов, получаются различные расчетные результаты.

function Start_interp1

x=[0.43 0.48 0.55 0.62 0.7 0.75];

y=[1.6359 1.73234 1.87686 2.03345 2.22846 2.35973];

n=length(x);

t=[0.45 0.52 0.59 0.67 0.72];

xi=[x(1):0.02:x(length(x))];

yiN=interp1(x,y,xi,'nearest');

yiL=interp1(x,y,xi,'linear');

yiS=interp1(x,y,xi,'spline');

yTS=interp1(x,y,t,'spline');

yTN=interp1(x,y,t,'nearest');

yTL=interp1(x,y,t,'linear');

plot(x,y,'o',xi,yiN,'r-',xi,yiL,'g-',xi,yiS,'c-',t,yTN,'x',t,yTL,'*',t,yTS,'+');

grid on;

end

linear линейная интерполяция;

spline кубический сплайн;

nearest интерполяция по соседним точкам (значение равно значению в ближайшей узловой точке).

Подбор параметров. Для подбора параметров предназначена функция [ res ]= lsqcurvefit (fun, a 0, x, y, Lx, Rx), где fun имя М-функции, вычисляющей левую часть выражения вида y = f (x, a ₀, a ₁, …, a_k), параметры которой a ₀, a ₁, …, a_k неизвестны и их-то и надо найти. Они перечислены в векторе возвращаемых значений res. а0 начальное приближение параметров. x, y экспериментальные данные. Lx, Rx ограничения на параметры вида Lx ≤ x ≤ Rx, Lx ≤ x, x ≤ Rx.

Пусть, для примера (Рис.1.37), надо подобрать параметры функции y = a ₁ + a ₂ x + a ₃ x ² + a ₄ x ³. Пусть даны экспериментальные данные.

Х	10.1	10.3	10.8	11	11.1	11.4	12.2	13.3	13.8	14	14.5	15	15.6	15.8	17	18.1	19
У	24	26	45	37	55	51	75	84	74	91	87	94	92	96	97	9 8	99

Выберем вектор начальных приближений а0=(2, 0, 1,0).

Опишем функцию exMnk, она вычисляет значение многочлена

function y=exMnk(a,x) y=a(1)+a(2)*x+a(3)*x.^2+a(4)*x.^3; end

Затем введем данные и построим график по ним.

>> X=[10.1,10.3,10.8,11,11.1,11.4,12.2,13.3,13.8,14,14.5,15,15.6,15.8,17,18.1,19]; >> Y=[24, 26, 45, 37, 55, 51, 75, 84, 74, 91, 87, 94, 92, 96, 97, 98, 99];

Применим функцию lsqcurvefit(@exMNK,a0,x,y), получим коэффициенты полинома.

>> plot(x,y,'o'),grid; >> a0=[2 0 1 0]; A=lsqcurvefit(@exMnk,a0,X,Y)		Maximum number of function evaluations exceeded; increase options.MaxFunEvals A = -19.6243 -8.5171 2.0875 -0.0689
function Start_lsq x=10:0.1:19; y=exMnk(A,x); plot(X,Y,'o'),grid; hold on; plot(x,y,'r-'); end	Зададим диапазон х по горизонтальной оси и построим (рис.1.37) график функции exMnk.

Рис.1.37. Экспериментальные данные и расчетный график.

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

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

Локальный минимум функции одной переменной ищет встроенная функция [x, y]=fminbnd(name, a, b, options). В ней:

· name имя М-функции,

· a, b границы поиска,

· options параметры управляющие ходом решения,

· x, y координаты точки минимума.

Для вычисления локального максимума, надо взять функцию name с другим знаком. Для примера, рассмотрим функцию y(x)=x⁴-0.5x³-28x²+140

Сначала найдем минимум на интервале, например, от -4 до -2.

[x,y]=fminbnd(@ext,-4,-2), где @ext обращение к М-функции вычисляющей интересующую нас функцию.

Потом найдем максимум на интервале от -6 до 6, для этого обратимся не к @ext, а к @ext_2, в которой функция y(x) взята с другим знаком, то есть y(x)=-(x⁴-0.5x³-28x²+140). То есть ее минимум будет максимумом исходной функции.

Получим и в том, и в другом случае координаты точки минимума (максимума) на заданном отрезке, что и требовалось.

function y=ext(x) y=x.^4-0.5*x.^3-28*x.^2+140; end function y=ext_2(x) y=-(x.^4-0.5*x.^3-28*x.^2+140); end function Start_fminb [x,y]=fminbnd(@ext,-6,6) [x,y]=fminbnd(@ext_2,-6,6) y=-y end

x = 3.9339 y = -84.2624 x = 5.0643e-006 y =  -140.0000 y = 140.0000