Метод Рунге — Кутты

Ме́тоды Ру́нге — Ку́тты (в литературе встречается название методы Рунге — Кутта) — большой класс численных методов решения задачи Коши для обыкновенных дифференциальных уравнений и их систем. Первые методы данного класса были предложены около 1900 года немецкими математиками К. Рунге и М. В. Куттой.

К классу методов Рунге — Кутты относятся явный метод Эйлера и модифицированный метод Эйлера с пересчётом, которые представляют собой соответственно методы первого и второго порядка точности. Существуют стандартные явные методы третьего порядка точности, не получившие широкого распространения. Наиболее часто используется и реализован в различных математических пакетах (Maple, MathCAD, Maxima) классический метод Рунге — Кутты, имеющий четвёртый порядок точности. При выполнении расчётов с повышенной точностью всё чаще применяются методы пятого и шестого порядков точности. Построение схем более высокого порядка сопряжено с большими вычислительными трудностями.

Методы седьмого порядка должны иметь по меньшей мере девять стадий, а методы восьмого порядка — не менее 11 стадий. Для методов девятого и более высоких порядков (не имеющих, впрочем, большой практической значимости) неизвестно, сколько стадий необходимо для достижения соответствующего порядка точности.

Классический метод Рунге — Кутты четвёртого порядкаПравить

Метод Рунге — Кутты четвёртого порядка при вычислениях с постоянным шагом интегрирования столь широко распространён, что его часто называют просто методом Рунге — Кутты.

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений первого порядка (далее $\mathbf{y},\mathbf{f},{\mathbf{k}}_i\in <!--\; \in \; -->{\mathbb{R}}^n$$\mathbf {y} ,\mathbf {f} ,\mathbf {k} _{i}\in \mathbb {R} ^{n}$ , а $x,h\in <!--\; \in \; -->{\mathbb{R}}^1$$x,h\in {\mathbb {R}}^{1}$ ).

${\mathbf{\text{y}}}^{\prime }=\mathbf{\text{f}}(x,\mathbf{\text{y}}),\mathbf{\text{y}}(x_0)={\mathbf{\text{y}}}_0.$${\textbf {y}}'={\textbf {f}}(x,{\textbf {y}}),\quad {\textbf {y}}(x_{0})={\textbf {y}}_{0}.$ 

Тогда приближенное значение в последующих точках вычисляется по итерационной формуле:

${\mathbf{\text{y}}}_{n+1}={\mathbf{\text{y}}}_n+\frac{h}{6}({\mathbf{\text{k}}}_1+2{\mathbf{\text{k}}}_2+2{\mathbf{\text{k}}}_3+{\mathbf{\text{k}}}_4)$$\textbf{y}_{n+1} = \textbf{y}_n + {h \over 6}(\textbf{k}_1 + 2\textbf{k}_2 + 2\textbf{k}_3 + \textbf{k}_4)$ 

Вычисление нового значения проходит в четыре стадии:

${\mathbf{\text{k}}}_1=\mathbf{\text{f}}\left(x_n,{\mathbf{\text{y}}}_n\right),$${\textbf {k}}_{1}={\textbf {f}}\left(x_{n},{\textbf {y}}_{n}\right),$ 

${\mathbf{\text{k}}}_2=\mathbf{\text{f}}\left(x_n+\frac{h}{2},{\mathbf{\text{y}}}_n+\frac{h}{2}{\mathbf{\text{k}}}_1\right),$${\textbf {k}}_{2}={\textbf {f}}\left(x_{n}+{h \over 2},{\textbf {y}}_{n}+{h \over 2}{\textbf {k}}_{1}\right),$ 

${\mathbf{\text{k}}}_3=\mathbf{\text{f}}\left(x_n+\frac{h}{2},{\mathbf{\text{y}}}_n+\frac{h}{2}{\mathbf{\text{k}}}_2\right),$${\textbf {k}}_{3}={\textbf {f}}\left(x_{n}+{h \over 2},{\textbf {y}}_{n}+{h \over 2}{\textbf {k}}_{2}\right),$ 

${\mathbf{\text{k}}}_4=\mathbf{\text{f}}\left(x_n+h,{\mathbf{\text{y}}}_n+h{\mathbf{\text{k}}}_3\right).$${\textbf {k}}_{4}={\textbf {f}}\left(x_{n}+h,{\textbf {y}}_{n}+h\ {\textbf {k}}_{3}\right).$ 

где $h$$h$  — величина шага сетки по $x$$x$ .

Этот метод имеет четвёртый порядок точности. Это значит, что ошибка на одном шаге имеет порядок $O(h^5)$$O(h^{5})$ , а суммарная ошибка на конечном интервале интегрирования имеет порядок $O(h^4)$$O(h^4)$  .

Явные методы Рунге — КуттыПравить

Семейство явных методов Рунге — Кутты является обобщением как явного метода Эйлера, так и классического метода Рунге — Кутты четвёртого порядка. Оно задаётся формулами

${\mathbf{\text{y}}}_{n+1}={\mathbf{\text{y}}}_n+h\underset{i=1}{\overset{n}{\sum <!--\; \sum \; -->}}{b}_i{\mathbf{\text{k}}}_i,$${\textbf {y}}_{n+1}={\textbf {y}}_{n}+h\sum _{i=1}^{n}b_{i}{\textbf {k}}_{i},$ 

где $h$$h$  — величина шага сетки по $x$$x$ , а вычисление нового значения проходит в $s$$s$  этапов:

${\begin{array}{ll}{\textbf {k}}_{1}=&{\textbf {f}}(x_{n},{\textbf {y}}_{n}),\\{\textbf {k}}_{2}=&{\textbf {f}}(x_{n}+c_{2}h,{\textbf {y}}_{n}+a_{21}h{\textbf {k}}_{1}),\\\cdots &\\{\textbf {k}}_{s}=&{\textbf {f}}(x_{n}+c_{s}h,{\textbf {y}}_{n}+a_{s1}h{\textbf {k}}_{1}+a_{s2}h{\textbf {k}}_{2}+\cdots +a_{s,s-1}h{\textbf {k}}_{s-1})\end{array}}$ 

Конкретный метод определяется числом $s$$s$  и коэффициентами $b_i,a_{ij}$$b_{i},a_{{ij}}$  и $c_i$$c_{i}$ . Эти коэффициенты часто упорядочивают в таблицу (называемую таблицей Бутчера):

${\begin{array}{c|ccccc}0&&&&&\\c_{2}&a_{{21}}&&&&\\c_{3}&a_{{31}}&a_{{32}}&&&\\\vdots &\vdots &\vdots &\ddots &&\\c_{s}&a_{{s1}}&a_{{s2}}&\dots &a_{{ss-1}}&\\\hline &b_{1}&b_{2}&\dots &b_{{s-1}}&b_{s}\end{array}}$ 

Для коэффициентов метода Рунге — Кутты должны быть выполнены условия $\underset{j=1}{\overset{i-<!--\; -\; -->1}{\sum <!--\; \sum \; -->}}{a}_{ij}=c_i$$\sum _{{j=1}}^{{i-1}}a_{{ij}}=c_{i}$  для $i=2,\dots <!--\; \dots \; -->,s$$i=2,\ldots ,s$ . Если требуется, чтобы метод имел порядок $p$$p$ , то следует также обеспечить условие

$\stackrel{¯<!--\; ¯\; -->}{\mathbf{\text{y}}}(h+x_0)-<!--\; -\; -->\mathbf{\text{y}}(h+x_0)=O(h^{p+1}),$${\bar {{\textbf {y}}}}(h+x_{0})-{{\textbf {y}}}(h+x_{0})=O(h^{{p+1}}),$ 

где $\stackrel{¯<!--\; ¯\; -->}{\mathbf{\text{y}}}(h+x_0)$${\bar {{\textbf {y}}}}(h+x_{0})$  — приближение, полученное по методу Рунге — Кутты. После многократного дифференцирования это условие преобразуется в систему полиномиальных уравнений относительно коэффициентов метода.

Неявные методы Рунге — КуттыПравить

Все до сих пор упомянутые методы Рунге — Кутты являются явными методами^[en]. К сожалению, явные методы Рунге — Кутты, как правило, непригодны для решения жестких уравнений из-за малой области их абсолютной устойчивости. Неустойчивость явных методов Рунге — Кутты создаёт весьма серьёзные проблемы и при численном решении дифференциальных уравнений в частных производных^[en].

Неустойчивость явных методов Рунге — Кутты мотивировала развитие неявных методов. Неявный метод Рунге — Кутты имеет вид:

$y_{n+1}=y_n+h\underset{i=1}{\overset{s}{\sum <!--\; \sum \; -->}}{b}_i{k}_i,$$y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i},$ 

где

$k_i=f(x_n+c_{i}h,y_n+h\underset{j=1}{\overset{s}{\sum <!--\; \sum \; -->}}{a}_{ij}{k}_j),i=1,\dots <!--\; \dots \; -->,s.$$k_{i}=f{\bigl (}x_{n}+c_{i}h,y_{n}+h\sum _{j=1}^{s}a_{ij}k_{j}{\bigr )},\quad i=1,\ldots ,s.$ 

Явный метод характерен тем, что матрица коэффициентов $a_{ij}$$a_{{ij}}$  для него имеет нижний треугольный вид (включая и нулевую главную диагональ) — в отличие от неявного метода, где матрица имеет произвольный вид. Это также видно по таблице Батчера^[en].

${\begin{array}{c|cccc}c_{1}&a_{{11}}&a_{{12}}&\dots &a_{{1s}}\\c_{2}&a_{{21}}&a_{{22}}&\dots &a_{{2s}}\\\vdots &\vdots &\vdots &\ddots &\vdots \\c_{s}&a_{{s1}}&a_{{s2}}&\dots &a_{{ss}}\\\hline &b_{1}&b_{2}&\dots &b_{s}\\\end{array}}={\begin{array}{c|c}{\mathbf {c}}&A\\\hline &{\mathbf {b^{T}}}\\\end{array}}$ 

Следствием этого различия является необходимость на каждом шагу решать систему уравнений для $k_i,i=1,2,...,s$$k_{i},i=1,2,...,s$ , где $s$$s$  — число стадий. Это увеличивает вычислительные затраты, однако при достаточно малом $h$$h$  можно применить принцип сжимающих отображений и решать данную систему методом простой итерации. В случае одной итерации это увеличивает вычислительные затраты всего лишь в два раза.

С другой стороны, Жан Кунцма́н^[fr] (1961) и Джон Батчер^[en] (1964) показали, что при любом количестве стадий $s$$s$  существует неявный метод Рунге — Кутты с порядком точности $p=2s$$p=2s$ . Это значит, например, что для описанного выше явного четырёхстадийного метода четвёртого порядка существует неявный аналог с вдвое большим порядком точности.

Неявный метод Рунге — Кутты второго порядкаПравить

Простейшим неявным методом Рунге — Кутты является модифицированный метод Эйлера «с пересчётом». Он задаётся формулой:

${\mathbf{y}}_{n+1}={\mathbf{y}}_n+h\frac{\mathbf{f}(x_n,{\mathbf{y}}_n)+\mathbf{f}(x_{n+1},{\mathbf{y}}_{n+1})}{2}$$\mathbf {y} _{n+1}=\mathbf {y} _{n}+h{\frac {\mathbf {f} (x_{n},\mathbf {y} _{n})+\mathbf {f} (x_{n+1},\mathbf {y} _{n+1})}{2}}$ .

Для его реализации на каждом шаге необходимы как минимум две итерации (и два вычисления функции).

Прогноз:

${\stackrel{~<!--\; ~\; -->}{\mathbf{y}}}_{n+1}={\mathbf{y}}_n+h\mathbf{f}(x_n,{\mathbf{y}}_n)$${\tilde {\mathbf {y} }}_{n+1}=\mathbf {y} _{n}+h\mathbf {f} (x_{n},\mathbf {y} _{n})$ .

Коррекция:

${\mathbf{y}}_{n+1}={\mathbf{y}}_n+h\frac{\mathbf{f}(x_n,{\mathbf{y}}_n)+\mathbf{f}(x_{n+1},{\stackrel{~<!--\; ~\; -->}{\mathbf{y}}}_{n+1})}{2}$$\mathbf {y} _{n+1}=\mathbf {y} _{n}+h{\frac {\mathbf {f} (x_{n},\mathbf {y} _{n})+\mathbf {f} (x_{n+1},{\tilde {\mathbf {y} }}_{n+1})}{2}}$ .

Вторая формула — это простая итерация решения системы уравнений относительно ${\mathbf{y}}_{n+1}$$\mathbf {y} _{n+1}$ , записанной в форме сжимающего отображения. Для повышения точности итерацию-коррекцию можно сделать несколько раз, подставляя ${\stackrel{~<!--\; ~\; -->}{\mathbf{y}}}_{n+1}={\mathbf{y}}_{n+1}$${\tilde {\mathbf {y} }}_{n+1}=\mathbf {y} _{n+1}$ . Модифицированный метод Эйлера «с пересчётом» имеет второй порядок точности.

УстойчивостьПравить

Преимуществом неявных методов Рунге — Кутты в сравнении с явными является их бо́льшая устойчивость, что особенно важно при решении жестких уравнений. Рассмотрим в качестве примера линейное уравнение y' = λy. Обычный метод Рунге — Кутты, применённый к этому уравнению, сведётся к итерации $y_{n+1}=r(h\mathrm{\lambda <!--\; \lambda \; -->})y_n$$y_{{n+1}}=r(h\lambda )\,y_{n}$ , с r, равным

$r(z)=1+z{b}^T(I-<!--\; -\; -->zA{)}^{-<!--\; -\; -->1}e=\frac{det(I-<!--\; -\; -->zA+ze{b}^T)}{det(I-<!--\; -\; -->zA)},$$r(z)=1+zb^{T}(I-zA)^{{-1}}e={\frac {\det(I-zA+zeb^{T})}{\det(I-zA)}},$ 

где $e$$e$  обозначает вектор-столбец из единиц. Функция $r$$r$  называется функцией устойчивости. Из формулы видно, что $r$$r$  является отношением двух полиномов степени $s$$s$ , если метод имеет $s$$s$  стадий. Явные методы имеют строго нижнюю треугольную матрицу $A,$$A,$  откуда следует, что $det(I-<!--\; -\; -->zA)=1,$$\det(I-zA)=1,$  и что функция устойчивости является многочленом.

Численное решение данного примера сходится к нулю при условии $|r(z)|<1$  с $z=h\mathrm{\lambda <!--\; \lambda \; -->}$$z=h\lambda$ . Множество таких $z$$z$  называется областью абсолютной устойчивости. В частности, метод называется A-устойчивым, если все $z$$z$  с ${\textrm {Re}}(z)<0$  находятся в области абсолютной устойчивости. Функция устойчивости явного метода Рунге — Кутты является многочленом, поэтому явные методы Рунге — Кутты в принципе не могут быть A-устойчивыми.

Если метод имеет порядок $p$$p$ , то функция устойчивости удовлетворяет условию $r(z)={\text{e}}^z+O(z^{p+1})$$r(z)={\textrm {e}}^{z}+O(z^{{p+1}})$  при $z\to <!--\; \to \; -->0$$z\to 0$ . Таким образом, представляет интерес отношение многочленов данной степени, приближающее экспоненциальную функцию наилучшим образом. Эти отношения известны как аппроксимации Паде. Аппроксимация Паде с числителем степени $m$$m$  и знаменателем степени $n$$n$  А-устойчива тогда и только тогда, когда $m\le <!--\; \le \; -->n\le <!--\; \le \; -->m+2.$$m\leq n\leq m+2.$ 

$s$$s$ -стадийный метод Гаусса — Лежандра имеет порядок $2s$$2s$ , поэтому его функция устойчивости является приближением Паде $m=n=s$$m=n=s$ . Отсюда следует, что метод является A-устойчивым. Это показывает, что A-устойчивые методы Рунге — Кутты могут иметь сколь угодно высокий порядок. В отличие от этого, порядок А-устойчивости методов Адамса не может превышать два.

ПроизношениеПравить

Согласно грамматическим нормам русского языка, фамилия Ку́тта склоняется, поэтому говорят: «Метод Ру́нге — Ку́тты». Правила русской грамматики предписывают склонять все фамилии (в том числе и мужские), оканчивающиеся на -а, -я, которым предшествует согласный. Единственное исключение — фамилии французского происхождения с ударением на последнем слоге типа Дюма́, Золя́. Однако иногда встречается несклоняемый вариант «Метод Ру́нге — Ку́тта» (например, в книге).

Пример решения на алгоритмических языках программированияПравить

$y^{″}+4y=\cos <!--\; \; -->3x,y(0)=0.8,y^{\prime }(0)=2,x\in <!--\; \in \; -->[0,1],h=0.1$$y''+4y=\cos 3x,\quad y(0)=0.8,\ y'(0)=2,\ x\in [0,1],\ h=0.1$

производя замену $y^{\prime }=z$$y'=z$  и перенося $4y$$4y$  в правую часть, получаем систему:

$\{\begin{array}{l}{y}^{\prime }=z=g(x,y,z)\\ z^{\prime }=\cos <!--\; \; -->(3x)-<!--\; -\; -->4y=f(x,y,z)\end{array}$${\begin{cases}y'=z=g(x,y,z)\\z'=\cos(3x)-4y=f(x,y,z)\end{cases}}$

код на Java для решения системы дифференциальных уравнений методом Рунге-Кутты

public class MainClass {

	public static void main(String[] args) {
		int k = 2;
		double Xo, Yo, Y1, Zo, Z1;
		double k1, k2, k4, k3, h;
		double q1, q2, q4, q3;
                /*
                 *Начальные условия
                 */
		Xo = 0;
		Yo = 0.8;
                Zo = 2;

		h = 0.1; // шаг

		System.out.println("\tX\t\tY\t\tZ");
		for(; r(Xo,2)<1.0; Xo += h){
			
			k1 = h * f(Xo, Yo, Zo);
			q1 = h * g(Xo, Yo, Zo);
			
			k2 = h * f(Xo + h/2.0, Yo + q1/2.0, Zo + k1/2.0);
			q2 = h * g(Xo + h/2.0, Yo + q1/2.0, Zo + k1/2.0);
			
			k3 = h * f(Xo + h/2.0, Yo + q2/2.0, Zo + k2/2.0);
			q3 = h * g(Xo + h/2.0, Yo + q2/2.0, Zo + k2/2.0);
			
			k4 = h * f(Xo + h, Yo + q3, Zo + k3);
			q4 = h * g(Xo + h, Yo + q3, Zo + k3);
			
			Z1 = Zo + (k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
			Y1 = Yo + (q1 + 2.0*q2 + 2.0*q3 + q4)/6.0;
			System.out.println("\t" + r(Xo + h, k) + "\t\t" + r(Y1 ,k) + "\t\t" + r(Z1 ,k));
			Yo = Y1;
			Zo = Z1;
		}
		
	}
        /**
         * функция для округления и отбрасывания "хвоста"
         */
	public static double r(double value, int k){
		return (double)Math.round((Math.pow(10, k)*value))/Math.pow(10, k);
	}
        /**
         * функции, которые получаются из системы
         */
	public static double f(double x, double y, double z){
		return (Math.cos(3*x) - 4*y);
	}
	public static double g(double x, double y, double z){
		return (z);
	}

}

Код на языке C#

using System;
using System.Collections.Generic;

namespace PRJ_RungeKutta
{
    /// <summary>
    /// Реализация метода Ру́нге — Ку́тты для обыкновенного дифференциального уравнения
    /// </summary>
    public abstract class RungeKutta
    {
        /// <summary>
        /// Текущее время
        /// </summary>
        public double t;  
        /// <summary>
        /// Искомое решение Y[0] — само решение, Y[i] — i-я производная решения
        /// </summary>
        public double[] Y; 
        /// <summary>
        /// Внутренние переменные 
        /// </summary>
        double[] YY, Y1, Y2, Y3, Y4;
        protected double[] FY;
        /// <summary>
        /// Конструктор
        /// </summary>
        /// <param name="N">размерность системы</param>
        public RungeKutta(uint N)  
        {
            Init(N);
        }
        /// <summary>
        /// Конструктор
        /// </summary>
        public RungeKutta(){}
        /// <summary>
        /// Выделение памяти под рабочие массивы
        /// </summary>
        /// <param name="N">Размерность массивов</param>
        protected void Init(uint N)
        {             
            Y = new double[N];
            YY = new double[N];
            Y1 = new double[N];
            Y2 = new double[N];
            Y3 = new double[N];
            Y4 = new double[N];
            FY = new double[N];
        }
        /// <summary>
        /// Установка начальных условий
        /// </summary>
        /// <param name="t0">Начальное время</param>
        /// <param name="Y0">Начальное условие</param>
        public void SetInit(double t0, double[] Y0) 
        {                                     
            t = t0;
            if (Y == null) 
                Init((uint)Y0.Length);
            for (int i = 0; i < Y.Length; i++)
                Y[i] = Y0[i];
        }
        /// <summary>
        /// Расчет правых частей системы
        /// </summary>
        /// <param name="t">текущее время</param>
        /// <param name="Y">вектор решения</param>
        /// <returns>правая часть</returns>
        abstract public double[] F(double t, double[] Y); 
        /// <summary>
        /// Следующий шаг метода Рунге-Кутта
        /// </summary>
        /// <param name="dt">текущий шаг по времени (может быть переменным)</param>
        public void NextStep(double dt)
        {
            int i;

            if (dt < 0) return;

            // рассчитать Y1
            Y1 = F(t, Y);

            for (i = 0; i < Y.Length; i++)
                YY[i] = Y[i] + Y1[i] * (dt / 2.0);
            
            // рассчитать Y2
            Y2 = F(t + dt / 2.0, YY);

            for (i = 0; i < Y.Length; i++)
                YY[i] = Y[i] + Y2[i] * (dt / 2.0);
            
            // рассчитать Y3
            Y3 = F(t + dt / 2.0, YY);

            for (i = 0; i < Y.Length; i++)
                YY[i] = Y[i] + Y3[i] * dt;
            
            // рассчитать Y4
            Y4 = F(t + dt, YY);

            // рассчитать решение на новом шаге
            for (i = 0; i < Y.Length; i++)
                Y[i] = Y[i] + dt / 6.0 * (Y1[i] + 2.0 * Y2[i] + 2.0 * Y3[i] + Y4[i]);

            // рассчитать текущее время
            t = t + dt; 
        }
    }
    class TMyRK : RungeKutta
    {
        public TMyRK(uint N) : base(N) { }

        /// <summary>
        /// пример математический маятник 
        /// y''(t)+y(t)=0
        /// </summary>
        /// <param name="t">Время</param>
        /// <param name="Y">Решение</param>
        /// <returns>Правая часть</returns>
        public override double[] F(double t, double[] Y)
        {
            FY[0] =  Y[1]; 
            FY[1] = -Y[0]; 
            return FY;
        }
        /// <summary>
        /// Пример использования
        /// </summary>
        static public void Test()
        {
            // Шаг по времени
            double dt = 0.001;
            // Объект метода
            TMyRK task = new TMyRK(2);
            // Определим начальные условия y(0)=0, y'(0)=1 задачи
            double[] Y0 = { 0, 1 }; 
            // Установим начальные условия задачи
            task.SetInit(0, Y0);
            // решаем до 15 секунд
            while (task.t <= 15) 
            {
                Console.WriteLine("Time = {0:F5}; Func = {1:F8};  d Func / d x = {2:F8}", task.t, task.Y[0], task.Y[1]); // вывести t, y, y'
                // рассчитать на следующем шаге, шаг интегрирования 
                task.NextStep(dt);
            }
            Console.ReadLine();
        }
    }

    class Program
    {
        static void Main(string[] args)
        {
            TMyRK.Test();
        }
    }
}

В программе на С# используется абстрактный класс RungeKutta, в котором следует переопределить абстрактный метод F, задающий правые части уравнений.

Пример решения в среде MATLABПравить

Решение систем дифференциальных уравнений методом Рунге — Кутты является одним из самых распространённых численных методов решений в технике. В среде MATLAB реализована его одна из разновидностей — метод Дормана — Принса^[en]. Для решения системы уравнений необходимо сначала записать функцию, вычисляющую производные, то есть функции y = g(x, y, z) и z = cos(3x) − 4y = f(x, y, z), о чём сказано выше. В одном из каталогов, к которому имеется доступ из системы MATLAB, нужно создать текстовый файл с именем (например) runge.m со следующим содержимым (для MATLAB версии 5.3):

MATLAB, runge.m

function Dy = runge(x, y)
Dy = y(:);
Dy(1) = y(2);
Dy(2) = cos(3*x) - 4*y(1);

Имя файла и имя функции должно совпадать, но оно может быть любым неиспользуемым ранее.

Затем необходимо создать главный файл c именем, например, main.m, который будет выполнять основные вычисления. Этот главный файл будет содержать следующий текст:

MATLAB, main.m

clear; clc; % Очистка памяти и экрана
h = 0.1; % Шаг интегрирования
x_fin = 8; % Конечное время интегрирования
y0 = 0.8; % Начальное значение функции
Dy0 = 2; % Начальное значение производной функции
[x, y] = ode45('runge', [0:h:x_fin], [y0 Dy0]); % Метод Рунге — Кутты
plot(x, y, 'LineWidth', 2); grid; % Построение графика и сетки
legend('y(x)', 'y''(x)', 0); % Легенда на графике

Так как MATLAB ориентирован на работу с матрицами, решение по методу Рунге — Кутты очень легко выполняется для целого ряда x, как, например, в приведённом примере программы. Здесь решение — график функции в пределах времён от 0 до x_fin.

 

Решение в среде MATLAB

Переменные x и y, полученные в результате работы функции ODE45, есть векторы значений. Очевидно, что решение конкретно заданного выше примера — второй элемент x, так как первое значение — 0, шаг интегрирования h = 0,1, а интересующее значение x = 0,1. Следующая запись в командном окне MATLAB даст искомое решение:

MATLAB, решение

y1 = y(find(x == 0.1))

Ответ: y1 = 0,98768

ЛитератураПравить

• Hairer E., Wanner G. . Solving ordinary differential equations II: Stiff and differential-algebraic problems. 2nd ed. — Berlin, New York: Springer-Verlag, 1996. — ISBN 978-3-540-60452-5.
• Iserles A. . A First Course in the Numerical Analysis of Differential Equations. — Cambridge: Cambridge University Press, 1996. — ISBN 978-0-521-55655-2.
• Süli E., Mayers D. . An Introduction to Numerical Analysis. — Cambridge: Cambridge University Press, 2003. — ISBN 0-521-00794-1.