Skip to content

Решение ОДУ в Matlab

Доброго времени суток! Сегодня мы поговорим о решении ОДУ (обыкновенных дифференциальных уравнений) в 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);

Погрешность метода 1 порядка
На данном графике показана зависимость величины ошибки от шага.

Метод Рунге-Кутта второго порядка

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».

Скачать исходник Рунге-Кутт первого порядка
Скачать исходник Рунге-Кутт второго порядка
Скачать исходник Рунге-Кутт четвёртого порядка

Опубликовано вMatlab

2 комментария

  1. Анна Анна

    Мне в курсовой попался метод Рунге-Кутта третьего порядка. Можете пожалуйста скинуть информацию о нем?

    • Brute Brute

      Здравствуй, Анна!

      Используй этот блок в цикле:
      x(n)=(n-1)*h;
      k1=f(x(n),y(n));
      k2=f(x(n)+h/3 ,y(n)+0.5*(h/3)*k1);
      k3=f(x(n)+(2*h)/3,y(n)+(2*h)/3*k2);
      y(n+1)=y(n)+(h/4)*(k1+3*k3)
      И следующую строку для проверки погрешности:
      M(k)=abs(y(N)-exp(0.5))/h^3;

      P.S. набросал без проверки, но должно работать. Отпишись, как проверишь 🙂

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *