Доброго времени суток! Сегодня мы поговорим о решении ОДУ (обыкновенных дифференциальных уравнений) в Matlab. Перед тем как мы начнём обсуждать данную тему, советую вам ознакомиться с темой: Численное дифференцирование в Matlab, чтобы лучше понимать теоретическую составляющую решения ОДУ.
Обыкновенные дифференциальные уравнения
С помощью дифференциальных уравнений можно описать разные задачи: движения системы, взаимодействующих материальных точек, химической кинетики и т.д. Различают три типа задач для систем диф. уравнений:
- Задача Коши
- Краевая задача
- Задача на собственные значения
Кратко расскажу о их сути:
Задача Коши предполагает дополнительные условия в виде значения функции в определённой точке.
Краевая задача подразумевает поиск решения на заданном отрезке с краевыми (граничными) условиями в концах интервала или на границе области.
Задача на собственные значения — помимо искомых функций и их производных, в уравнение входят дополнительное несколько неизвестных параметров, которые являются собственными значениями.
Методы решения дифференциальных уравнений
Решение ОДУ в Matlab и не только, в первую очередь, сводится к выбору порядка численного метода решения. Порядок численного метода не связан с порядком дифференциального уравнения. Высокий порядок у численного метода означает его скорость сходимости.
В случае большого интервала, с помощью алгоритмов с низким порядком сжимают интервал с решениями и находят приблизительные корни, а затем уже уточняют корни с помощью методов с высоким порядком.
Решение обыкновенных дифференциальных уравнений в Matlab можно реализовать «своими ручками», прописав алгоритм по разным схемам. Но также в Matlab есть встроенные функции, выполняющие все стандартные задачи.
Метод Рунге-Кутта первого порядка
yi+1 = yi + h*f(xi, yi)
Методы Рунге-Кутта представляют собой разложения в ряд Тейлора и от количества использованных элементов ряда зависит порядок этого метода. Следовательно, помимо Рунге-Кутта первого порядка, вы сможете увидеть методы других порядков. Иногда их называют другими именами.
Например, Метод Рунге-Кутта первого порядка, также известен как Метод Эйлера или Метод ломаных. Информацию о его математическом и графическом представлении советую поискать в гугл. Мы же поговорим о том, как Метод Рунге-Кутта первого порядка реализуется в Matlab для решения ОДУ. Например:
Решить и привести график ошибки уравнения y' = y*x методом Рунге-Кутта первого порядка (Методом Эйлера, Методом ломаных).
Setka=10:10:10000; for k=1:length(Setka) %определяем параметры сетки N=Setka(k); h=1.0/(N-1); %задаем начальное условие y(1)=1; %применяем алгоритм метода ломаных/Эйлера for n=1:(N-1) y(n+1)=y(n)+ h*((n-1)*h*y(n)); %где (n-1)*h -> x end %вычисляем величину M(1) в оценке %погрешности численного решения M(k)=abs(y(N)-exp(0.5))/h; step(k)=h; end %рисуем зависимость величины M(1) от шага plot(step,M);
На данном графике показана зависимость величины ошибки от шага.
Метод Рунге-Кутта второго порядка
yi+1 = yi + (h/2)*f(xi, yi) + (h/2)*f(xi+1, yi+1)
Также известен как Метод Эйлера-Коши. Как видите, во второй части уравнения происходит обращения к следующему шагу. Но как тогда быть, если нам ещё не известен следующий шаг? Всё просто. Метод Рунге-Кутта второго порядка — это всё тот же метод первого порядка, однако, на половине шага происходит нахождение «первичного» решения, а затем происходит его уточнение. Это позволяет поднять порядок скорости сходимости до двух.
Решить и привести график ошибки уравнения u' = u*x методом Рунге-Кутта второго порядка.
f=@(x,u)x*u; %задаем набор сеток Setka=10:50:10000; %организуем цикл расчетов с разными сетками for k=1:length(Setka) N=Setka(k); h=1.0/(N-1); %задаем начальное условие y(1)=1; %вычисляем приближенные значения решения %дифференциального уравнения for n=1:(N-1) x(n)=(n-1)*h; y(n+1)=y(n)+h*(0.75*f(x(n),y(n))+... 0.25*f(x(n)+h/0.5,y(n)+... (h/0.5)*f(x(n),y(n)))); end %находим константу M в оценке погрешности %приближенного решения |y(N)-u(N)|<=Mh^2, %она не должна зависеть от h M(k)=abs(y(N)-exp(0.5))/h^2; step(k)=h; end %рисуем зависимость константы M от шага сетки h plot(step,M);
По сравнению с Рунге-Куттом первого порядка изначальная ошибка уже гораздо меньше.
Мы не будем говорить о третьем порядке, потому что задачи на третий порядок встречаются редко, но если будет необходимо, пишите в комментариях, выложу.
Метод Рунге-Кутта четвёртого порядка
T1 = h*f(xi, yi)
T2 = h*f(xi + h/2, yi + T1/2)
T3 = h*f(xi + h/2, yi + T2/2)
T4 = h*f(xi + h/2, yi + T3)
yi+1 = yi + (T1 + 2*T2 + 2*T3 + T4)/6
Метод Рунге-Кутта четвёртого порядка считается самым распространённым. Тем не менее, работает он аналогично второму и третьему порядку.
Решить и привести график ошибки уравнения u' = u*x методом Рунге-Кутта четвёртого порядка.
f=@(x,u)x*u; %задаем набор сеток Setka=10:50:10000; %организуем цикл расчетов с разными сетками for k=1:length(Setka) N=Setka(k); h=1.0/(N-1); %задаем начальное условие y(1)=1; %вычисляем приближенные решения %дифференциального уравнения for n=1:(N-1) x(n)=(n-1)*h; k1=f(x(n),y(n)); k2=f(x(n)+0.5*h,y(n)+0.5*h*k1); k3=f(x(n)+0.5*h,y(n)+0.5*h*k2); k4=f(x(n)+h,y(n)+h*k3); y(n+1)=y(n)+(h/6)*(k1+2*k2+2*k3+k4); end %находим константу M в оценке погрешности %приближенного решения |y(N)-u(x(N))|<=Mh^4, %она не должна зависеть от h M(k)=abs(y(N)-exp(0.5))/h^4; step(k)=h; end %рисуем зависимость константы M от шага сетки h loglog(step,M);
Как видите, на последней картинке размерность ошибки на столько мала, что пришлось воспользоваться loglog() для лучшей видимости.
Решение ОДУ в Matlab стандартными средствами
Стоит отметить, что мы с вами разобрали только один самый известный метод решения ОДУ с разными порядками. Однако, методов очень много.
Для решения дифференциальных уравнений и систем в MATLAB предусмотрены следующие функции:
ode45 (f, interval, X0, [options])
ode23 (f, interval, X0, [options])
ode113 (f, interval, X0, [options])
ode15s (f, interval, X0, [options])
ode23s (f, interval, X0, [options])
ode23t (f, interval, X0, [options])
ode23tb (f, interval, X0, [options])
Входными параметрами этих функций являются:
- f — вектор-функция для вычисления правой части уравнения системы уравнений;
- interval — массив из двух чисел, определяющий интервал интегрирования дифференциального уравнения или системы;
- Х0 — вектор начальных условий системы дифференциальных уравнений;
- options — параметры управления ходом решения дифференциального уравнения или системы.
Все функции возвращают:
- массив Т — координаты узлов сетки, в которых ищется решение;
- матрицу X, i-й столбец которой является значением вектор-функции решения в узле Тi.
В функции ode45 реализован метод Рунге-Кутта 4-5 порядка точности, в функции ode23 также реализован метод Рунге-Кутта, но 2-3 порядка, а функция ode113 реализует метод Адамса.
Для решения жёстких систем предназначены функция ode15s, в которой реализован метод Гира, и функция ode23s, реализующая метод Розенброка. Для получения более точного решения жёсткой системы лучше использовать функцию ode15s. Для решения системы с небольшим числом жёсткости можно использовать функцию ode23t, а для грубой оценки подобных систем служит функция ode23tb.
Символьное решение обыкновенных дифференциальных уравнений произвольного порядка осуществляет функция dsolve r = dsolve(‘eq1,eq2,…’, ‘cond1,cond2,…‘, ‘v’)
Пример использования:
% в отдельном М-файле с именем pr10 function f=pr10(x,y) f=sin(x+y)+(3/2)*(x-y); end) %в основном М-файле ode113(@pr10,[0 20],0) %Метод Адамса: @pr10 – ссылка на М-функцию, % [0 20]- интервал, 0 - условие: y(0)=0
%пример (u'=u^2) expl=dsolve('Du=u^2','x')
На этом мы закончим. Если остались вопросы, задавайте их в комментариях. Также вы можете скачать исходники чтобы лучше понять тему: «Решение ОДУ в Matlab».
Скачать исходник Рунге-Кутт первого порядка
Скачать исходник Рунге-Кутт второго порядка
Скачать исходник Рунге-Кутт четвёртого порядка
Будьте первым, кто оставит комментарий