Компьютерное моделирование и решение нелинейных уравнений

Метод прямоугольников

Словесный алгоритм метода прямоугольников:

Весь участок [a,b] делим на n равных частей с шагом h=(b-a)/n.

Определяем значение yi подынтегральной функции f(x) в каждой части деления, т.е.$$y_i=f(x_i), i=\overline{0,n}.$$

В каждой части деления подынтегральную функцию f(x) аппроксимируем интерполяционным многочленом степени n = 0, т.е. прямой, параллельной оси OX. В результате вся подынтегральная функция на участке [a,b] аппроксимируется ломаной линией.

Для каждой части деления определяем площадь Si частичного прямоугольника.

Суммируем эти площади. Приближенное значение интеграла I равно сумме площадей частичных прямоугольников.

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

$$I=\sum \limits_{i=0}^{n-1}S_i = \sum \limits_{i=0}^{n-1}h \cdot y_i = h \cdot \sum \limits_{i=0}^{n-1}y_i.$$ (рис 12.3) Метод левых прямоугольников

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

$$I=\sum \limits_{i=1}^{n}S_i = \sum \limits_{i=1}^{n}h \cdot y_i = h \cdot \sum \limits_{i=1}^{n}y_i.$$ (рис 12.4) Метод правых прямоугольников

Точность каждого метода прямоугольников имеет порядок h.

Алгоритм вычисления интеграла построим в виде итерационного процесса поиска с автоматическим выбором шага. На каждом шаге будем уменьшать шаг в два раза, то есть увеличивать число шагов n в два раза. Выход из процесса поиска организуем по точности вычисления интеграла. Начальное число шагов n=2.Схема алгоритма методов прямоугольников представлена на рис.12.5.

(рис 12.5) Схема алгоритма метода прямоугольников (с автоматическим выбором шага)

Условные обозначения:

a,b - концы интервала,

$$\varepsilon$$ - заданная точность,

с=0 - метод левых прямоугольников,

с=1 - метод правых прямоугольников,

S1 - значение интеграла на предыдущем шаге,

S - значение интеграла на текущем шаге.

Метод трапеций

Словесный алгоритм метода трапеций:

Интервал [a,b] делим на n равных частей с шагом h=(b-a)/n.

Вычисляем значение подынтегральной функции в каждой узловой точке$$y_i=f(x), i=\overline{0,n}.$$

На каждом шаге подынтегральную функцию f(x) аппроксимируем прямой, соединяющей две соседние узловые точки. В результате вся подынтегральная функция на участке [a,b] заменяется ломаной линией проходящей через все узловые точки.

Вычисляем площадь каждой частичной трапеции.

Приближенное значение интеграла равно сумме площадей частичных трапеций, т.е.$$I=\sum \limits_{i=1}^{n-1}S_i.$$

Найдем площади Si частичных трапеций:

$$S_0=\frac{1}{2}h(y_0+y_1),\\ S_1=\frac{1}{2}h(y_1+y_2),\\ S_2=\frac{1}{2}h(y_2+y_3),\\ \ldots\\ S_{n-1}=\frac{1}{2}h(y_{n-1}+y_n).$$

Приближенное значение интеграла равно

$$I=\sum \limits_{i=1}^{n-1}S_i = \frac{h}{2} \sum \limits_{i=1}^{n}(y_i + y_{i+1}) = \frac{h}{2}(y_0 + y_n + 2\sum \limits_{i=1}^{n-1}y_i).$$

Точность метода трапеций имеет порядок h2.

Схема алгоритма метода трапеций представлена на рис.12.6.

(рис 12.6) Схема алгоритма метода трапеций (с автоматическим выбором шага)

Метод Симпсона

В методе Симпсона в каждой части деления подынтегральная функция аппроксимируется квадратичной параболой a0x2+a1x+a2. В результате вся кривая подынтегральной функции на участке [a,b] заменяется кусочно-непрерывной линией, состоящей из отрезков квадратичных парабол. Приближенное значение интеграла I равно сумме площадей под квадратичными параболами.

Т.к. для построения квадратичной параболы необходимо иметь три точки, то каждая часть деления в методе Симпсона включает два шага, т.е.

Lk=2h.

В результате количество частей деления N2=n/2. Тогда n в методе Симпсона всегда четное число.

Определим площадь S1 на участке [x0, x2] (рис.12.2).

Исходя из геометрического смысла определенного интеграла, площадь S1 равна определенному интегралу от квадратичной параболы на участке [x0, x2]:

$$S_1=\int_{x0}^{x2}(a_0x^2 + a_1x + a_2)dx = \frac{1}{3}a_0x^3 + \frac{1}{2}a_1x^2 + a_2x\left|_{xo}^{x2} \right. =\\ = \frac{x_2-x_0}{6}(2a_0(x_0 + x_0x_2 +x_2^2) + 3a_1(x_0 + x_2 + 6a^2).$$

Неизвестные коэффициенты квадратичной параболы а0 , а1, а2 определяем из условия прохождения параболой через три узловых точки с координатами (x0y0), (x1y1), (x2y2).

На основании этого условия строим систему линейных уравнений:

$$\left\{ \begin{array}{l} a_0x_0^2 + a_1x_0 + a_2 = y_0,\\ a_0x_1^2 + a_1x_1 + a_2 = y_1,\\ a_0x_2^2 + a_1x_2 + a_2 = y_2. \end{array} \right.$$

Решая эту систему, найдем коэффициенты параболы.

В результате имеем: $$S_1=\frac{h}{3}(y_0+4y_1+y_2)$$..

Для участка [x2, x4]: $$S_2=\frac{h}{3}(y_2+4y_3+y_4)$$..

:::::::::::::::::::

Для участка [xi-1, xi+1]: $$S_k=\frac{h}{3}(y_{i-1}+4y_i+y_{i+1})$$.,

где $$k=\frac{i+1}{2}$$.

Суммируя все площади S1 под квадратичными параболами, получим квадратурную формулу по методу Симпсона:

$$S=\sum \limits_{k=1}^{N2}S_k = \frac{h}{3} \sum \limits_{k=1}^{N2}(y_{i-1} + 4y_i + y_{i+1}),$$

где

N2 - количество частей деления.

Точность метода Симпсона имеет порядок (h3/h4).

Схема алгоритма метода Симпсона представлена на рис.12.7.

(рис 12.7) Схема алгоритма Симпсона (с автоматическим выбором шага)

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

Общий вид дифференциального уравнения

$$F(x,y,y')=0.$$

Нормальная форма дифференциального уравнения

$$y'=f(x,y),$$

где

y=y(x) -неизвестная функция, подлежащая определению,

f(x,y) - правая часть дифференциального уравнения в нормальной форме, равная первой производной функции y(x). В функцию f(x,y) помимо аргумента x входит и сама неизвестная функция y(x).

Пример:

$$x \cdot y'-(x^2-1) \cdot y=0$$ - общий вид дифференциального уравнения первого порядка,

$$y'=\frac{x^2-1}{x} \cdot y$$ - нормальная форма этого же уравнения.

Если неизвестная функция у зависит от одного аргумента x, то дифференциальное уравнение вида$$y'=f(x,y),$$ называется обыкновенным дифференциальным уравнением.

Если функция у зависит от нескольких аргументов, то такое дифференциальное уравнение называется дифференциальным уравнением в частных производных.

Общим решением обыкновенного дифференциального уравнения$$y'=f(x,y),$$ является семейство функций у=у(х,с) (рис 12.8):

(рис 12.8)

При решении прикладных задач ищут частные решения дифференциальных уравнений. Выделение частного решения из семейства общих решений осуществляется с помощью задания начальных условий:$$y\left|_{x=x_0} = y_0 \right.,$$ т.е. начальной точки с координатами (х₀, у₀).

Нахождение частного решения дифференциального уравнения$$y'=f(x,y),$$ удовлетворяющего начальному условию$$y\left|_{x=x_0} = y_0 \right.,$$ называется задачей Коши.

В численных методах задача Коши ставится следующим образом: найти табличную функцию $$y_i=f(x_i), i=\overline{1,n}$$ которая удовлетворяет заданному дифференциальному уравнению (12.2) и начальному условию (12.3) на отрезке [a,b] с шагом h, то есть найти таблицу

i	x	y
0	x0	y0
1	x1	y1
2	x2	y2
3	x3	y3
...	...	...
n	xn	yn

Здесь

h - шаг интегрирования дифференциального уравнения,

a=x0 - начало участка интегрирования уравнения,

b=xn - конец участка,

n=(b-a)/h - число шагов интегрирования уравнения.

На графике (рис 12.9) решение задачи Коши численными методами представляется в виде совокупности узловых точек с координатами (xi ,yi), $$i=\overline{1,n}$$.

(рис 12.9)

Методы Рунге - Кутта

Наиболее эффективными и часто встречаемыми методами решениями задачи Коши являются методы Рунге - Кутта. Они основаны на аппроксимации искомой функции у(х) в пределах каждого шага многочленом, который получен при помощи разложения функции у(х) в окрестности шага h каждой i-ой точки в ряд Тейлора:

$$y(x_i+h) = y(x_i) + h \cdot y'(x_i) + \frac{h^2}{2!}y''(x_i) + \\ + \frac{h^3}{3!}y'''(x_i) + \frac{h^4}{4!}y^{(4)}(x_i) + \frac{h^5}{5!}y^{(5)}(x_i) + \ldots$$

Усекая ряд Тейлора в различных точках и отбрасывая правые члены ряда, Рунге и Кутта получали различные методы для определения значений функции у(х) в каждой узловой точке. Точность каждого метода определяется отброшенными членами ряда.

Метод Рунге - Кутта 1-го порядка (метод Эйлера)

Отбросим в (12.4) члены ряда, содержащие h2, h3, h4.

Тогда $$y(x_i+h)=y(x_i)+h \cdot y'(x_i)$$.

Так как $$y'(x_i)=f(x_i,y_i)$$

Получим формулу Эйлера:

$$y_{i+1}=y_i + h \cdot f(x_i,y_i)$$

Так как точность методов Рунге-Кутта определяется отброшенными членами ряда (12.4), то точность метода Эйлера на каждом шаге составляет $$\approx h^2$$.

Алгоритм метода Эйлера можно построить в виде двух программных модулей: основной программы и подпрограммы ELER, реализующей метод.

(рис 12.10) Схема алгоритма метода Эйлера

Здесь

(x,y) -при вводе начальная точка, далее текущие значения табличной функции,

h -шаг интегрирования дифференциального уравнения,

b -конец интервала интегрирования.

Рассмотрим геометрический смысл метода Эйлера.

Формула Эйлера имеет вид:

$$y_{i+1}=y_i + h \cdot f(x_i,y_i),$$

где $$f(x_i,y_i) = y'(x_i) = \tg\alpha_i$$.

Тогда формула Эйлера принимает вид:

$$y_{i+1}=y_i + h \cdot \tg\alpha_i,$$

где

$$\tg\alpha_i$$ - тангенс угла наклона касательной к искомой функции у(x) в начальной точке каждого шага.

(рис 12.11) Геометрический смысл метода Эйлера

В результате в методе Эйлера на графике (рис 12.10) вся искомая функция y(x) на участке [a,b] аппроксимируется ломаной линией, каждый отрезок которой на шаге h линейно аппроксимирует искомую функцию. Поэтому метод Эйлера получил еще название метода ломаных.

В методе Эйлера наклон касательной в пределах каждого шага считается постоянным и равным значению производной в начальной точке шага xi. В действительности производная, а, значит, и тангенс угла наклона касательной к кривой y(x) в пределах каждого шага меняется. Поэтому в точке xi+h наклон касательной не должен быть равен наклону в точке xi. Следовательно, на каждом шаге вносится погрешность.

Первый отрезок ломаной действительно касается искомой интегральной кривой y(x) в точке (x0,y0). На последовательных же шагах касательные проводятся из точек (xi,yi), подсчитанных с погрешностью. В результате с каждым шагом ошибки накапливаются.

Основной недостаток метода Эйлера - систематическое накопление ошибок. Поэтому метод Эйлера рекомендуется применять для решения дифференциальных уравнений при малых значениях шага интегрирования h.

Метод Рунге - Кутта 2-го порядка (модифицированный метод Эйлера)

Отбросим в (12.4) члены ряда, содержащие h3, h4, h5:.

Тогда

$$y(x_i+h)=y(x_i)+h\cdot y'(x_i)+\frac{h^2}{2!} \cdot y"(x_i).$$

Чтобы сохранить член ряда, содержащий h2, надо определить вторую производную y"(xi).Ее можно аппроксимировать разделенной разностью 2-го порядка

$$y"=(x_i)=\frac{\Delta y'}{\Delta x}= \frac{y'(x_i+h)-y'(x_i)}{h}$$

Подставляя это выражение в (12.6), получим

$$y=(x_i+h) =y(x_i) + h\cdot y'(x_i) + \frac{h}{2} \cdot \frac{y'(x_i+h)-y'(x_i)}{h}=\\ =y(x_i) + \frac{h}{2}y'(x_i) + \frac{h}{2}y'(x_i+h).$$

Окончательно, модифицированная или уточненная формула Эйлера имеет вид:

$$y_{i+1}=y_i + \frac{h}{2}f(x_i,y_i) + \frac{h}{2}f(x_{i+1},y_{i+1}).$$

Как видно, для определения функции y(x) в точке i+1 необходимо знать значение правой части дифференциального уравнения f(xi+1, yi+1) в этой точке, для определения которой необходимо знать предварительное значение yi+1.

Для определения предварительного значения yi+1 воспользуемся формулой Эйлера. Тогда все вычисления на каждом шаге по модифицированной или уточненной формуле Эйлера будем выполнять в два этапа:

На первом этапе вычисляем предварительное значение $$y_{i+1}^Э$$ по формуле Эйлера

$$y_{i+1}^Э = y_i + hf(x_i,y_i).$$

На втором этапе уточняем значение y=i+1 по модифицированной или уточненной формуле Эйлера

$$y_{i+1} = y_i + \frac{h}{2}f(x_i,y_i) + \frac{h}{2}f(x_{i+1},y_{i+1}^Э).$$

Точность метода определяется отброшенными членами ряда Тейлора (12.4), т.е. точность уточненного или модифицированного метода Эйлера на каждом шаге $$\approx h^3$$.

Рассмотрим геометрический смысл модифицированного метода Эйлера.

Так как$$f(x_i,y_i) = y'(x) = \tg\alpha_i,\\ f(x_{i+1},y_{i+1}) = y'(x_{i+1}) = \tg\alpha_{i+1},$$ то модифицированную формулу Эйлера можно представить в виде:

$$y_{i+1} = y_I + \frac{h}{2}\tg\alpha_i + \frac{h}{2}\tg\alpha_{i+1},$$

где

$$\tg \alpha_i$$ - тангенс угла наклона касательной к искомой функции у(х) в начальной точке каждого шага,

$$\tg \alpha_{i+1}$$ - тангенс угла наклона касательной к искомой функции у(х) в конечной точке каждого шага.

(рис 12.12) Геометрический смысл модифицированного метода Эйлера

Здесь

P1 - накопленная ошибка в (i+1)й точке по методу Эйлера,

P2 - накопленная ошибка в (i+1)й точке по модифицированному методу Эйлера.

Как видно из рис.12.11, в первой половине каждого шага, то есть на участке [xi, xi+h/2], искомая функция y(x) аппроксимируется прямой, которая выходит из точки (xi, yi) под углом, тангенс которого $$\tg \alpha_i = f(x_i,y_i)$$.

Во второй половине этого же шага, т.е. на участке [xi + h/2,xi + h], искомая функция y(x) аппроксимируется прямой, которая выходит из точки с координатами

$$x=x_i + h/2,\\ y=y_i = \frac{h}{2} \cdot f(x_i,y_i)$$

под углом, тангенс которого

$$\tg \alpha_{i+1} = f(x_{i+1},y_{i+1}^Э).$$

В результате в модифицированном методе Эйлера функция у(х) на каждом шаге аппроксимируется не одной прямой, а двумя.

Алгоритм модифицированного метода Эйлера можно построить в виде двух программных модулей: основной программы и подпрограммы МELER, реализующей метод (рис. 12.13).

(рис 12.13) Схема алгоритма модифицированного метода Эйлера

Здесь

(x,y) -при вводе начальная точка, далее текущие значения табличной функции,

h -шаг интегрирования дифференциального уравнения,

b -конец интервала интегрирования.

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

Самое большое распространение из всех численных методов решения дифференциальных уравнений с помощью ЭВМ получил метод Рунге-Кутта 4-го порядка. В литературе он известен как метод Рунге-Кутта.

В этом методе на каждом шаге интегрирования дифференциальных уравнений искомая функция y(x) аппроксимируется рядом Тейлора (12.4), содержащим члены ряда с h4:

$$y(x_i+h) = y(x_i) + h \cdot y'(x_i) + \frac{h^2}{2!}y''(x_i) + \frac{h^3}{3!}y'''(x_i) + \frac{h^4}{4!}y^{(4)}(x_i)\ldots$$

В результате ошибка на каждом шаге имеет порядок h5.

Для сохранения членов ряда, содержащих h2,h3,h4 необходимо определить вторую y", третью y"' и четвертую y(4) производные функции y(x). Эти производные аппроксимируем разделенными разностями второго, третьего и четвертого порядков соответственно.

В результате для получения значения функции yi+1 по методу Рунге-Кутта выполняется следующая последовательность вычислительных операций:

$$T_1=h \cdot f(x_i,y_i),\\ T_2=h \cdot f(x_i+h/2,y_i+T_1/2),\\ T_3=h \cdot f(x_i+h/2,y_i+T_2/2),\\ T_4=h \cdot f(x_i+h/2,y_i+T_3),\\ y_{i+1}=y_1+(T_1+2 \cdot T_2+2 \cdot T_3+T_4)/6.$$

Вывод формулы не приведен. Предоставляется возможность вывод формул выполнить самостоятельно.

Алгоритм метода Рунге-Кутта (4-го порядка) можно построить в виде двух программных модулей: основной программы и подпрограммы Rk4, реализующей метод (рис 12.14).

(рис 12.14) Схема алгоритма метода Рунге-Кутта 4-го порядка.

Здесь

(x,y) -при вводе начальная точка, далее текущие значения табличной функции,

h -шаг интегрирования дифференциального уравнения,

b -конец интервала интегрирования.

Решение дифференциальных уравнений высоких порядков

Методы Рунге-Кутта можно использовать не только для решения дифференциальных уравнений первого порядка$$y'=f(x,y)$$ но и для решения дифференциальных уравнений более высоких порядков

$$y''=f(x,y,y'),\\ y'''=f(x,y,y',y''),\\ \ldots\\ y^{(m)}=f(x,y,y',y''',\ldots,y^{(m-1)}).$$

Любое дифференциальное уравнение m-го порядка

$$y^{(m)}=f(x,y,y',y",\ldots,y^{(m-1)})$$

можно свести к системе, состоящей из m уравнений первого порядка при помощи замен.

Заменим:

$$y_1=y',\\ y_2=y''=y'_1,\\ y_3=y'''=y'_2,\\ \ldots\\ y_m=y^{(m)}=y'_{(m-1)}.$$

В результате дифференциальное уравнение m -го порядка (12.8) сводится к системе, состоящей из m дифференциальных уравнений первого порядка:

$$\left\{ \begin{array}{l} y'=y_1,\\ y'_1=y_2,\\ y'_2=y_3,\\ \ldots\\ y'_{m-1}=f(x,y,y_1,y_2,\ldots,y_{m-1}). \end{array} \right.$$

Решением системы (12.2), а значит и дифференциального уравнения m -го порядка (12.1) является m табличных функций $$y, y_1=y', y_2=y_1^{''}, \ldots, y_m=y_{(m-1)}$$

Решение дифференциальных уравнений второго порядка

В задачах моделирования динамических систем наиболее часто приходится решать дифференциальные уравнения второго порядка.

Общий вид дифференциальных уравнений второго порядка:

$$F(x,y,y',y'')=0$$

Нормальная форма дифференциальных уравнений второго порядка:

$$y''=f(x,y,y')$$

Пример

Уравнение в общем виде

$$x^2y'' - xy' + (x^2 - 1)y = 0.$$

Его нормальная форма

$$y''=\frac{1}{x}y' - \frac{x^2-1}{x^2}y.$$

Дифференциальное уравнение второго порядка (12.11) можно свести к системе, состоящей из двух дифференциальных уравнений первого порядка при помощи замен.

Заменим y1=y',

Тогда y'1=y".

В результате уравнение (12.11) сводится к системе, состоящей из двух дифференциальных уравнений первого порядка:

$$\left\{ \begin{array}{l} y_1=y'\\ y'_1=f(x,y,y_1). \end{array} \right.$$

Для примера (12.12) эта система имеет вид:

$$\left\{ \begin{array}{l} y'=y_1\\ y_1^{'}=\frac{1}{x}y_1-\frac{x^2-1}{x^2}y \end{array} \right.$$

Решением этой системы являются две функции y(x) и y1(x),

где

$$y_1(x)=y'(x).$$

Сформулируем задачу Коши для системы, состоящей из двух дифференциальных уравнений второго порядка.

Дана система

$$\left\{ \begin{array}{l} y'=f_1(x,y,y_1),\\ y_1^{'}=f_2(x,y,y_1). \end{array} \right.$$

Даны два начальных условия:

$$y(x_0) = y_0,\\ y_1(x_0)=(y_1)_0.$$

Необходимо проинтегрировать систему на участке [a, b] с шагом h.

В численных методах задача Коши для системы (12.14) ставится следующим образом:

Найти табличные функции

$$y_i(x_i)$$ и $$(y_1)_i(x_i), i=\overline{1,n}$$, т.е. найти таблицу

i	x	y	y1
0	x0	y0	(y1)0
1	x1	y1	(y1)1
2	x2	y2	(y1)2
3	x3	y3	(y1)3
...	...	...	...
n	xn	yn	(y1)n

Здесь

h - шаг интегрирования дифференциального уравнения,

a=x0 - начало участка интегрирования уравнения,

b=xn - конец участка,

n=(b-a)/h - число шагов интегрирования уравнения.

На графике решением задачи Коши для системы, состоящей из двух дифференциальных уравнений первого порядка, является совокупность узловых точек (рис. 12.15).

При этом на каждом шаге, т.е. для каждого значения xi решением являются две узловые точки с координатами (xi, yi), (xi, (y1)i).

(рис 12.15)

Для решения системы дифференциальных уравнений используем те же методы, что и для решения одного дифференциального уравнения первого порядка. При этом необходимо соблюдать условие: на каждом шаге интегрирования, т. е. в точках с координатами х1 , х2 , х3 , : , хn все уравнения системы надо решать параллельно.

Для вычисления правых частей уравнений системы (12.14) необходимо сформировать подпрограмму PRAV.

Вернемся к примеру (12.13). Здесь на каждом шаге в подпрограмме PRAV будем вычислять правые части каждого уравнения системы:

Схема алгоритма решения системы (12.13) представлена на рис 12.16.

(рис 12.16) Схема алгоритма решения системы (12.6)

Здесь

h - шаг интегрирования дифференциального уравнения,

b - конец участка,

n - число шагов интегрирования уравнения,

x, y, y1 - при вводе начальные значения, далее - текущие значения табличной функции.

Решение дифференциальных уравнений m-го порядка методом Рунге-Кутта (4-го порядка)

Как уже было сказано, любое дифференциальное уравнение m-го порядка

$$y_{(m)}=f(x,y,y',y'',\ldots,y^{(m-1)})$$

сводится к системе, состоящей из m дифференциальных уравнений 1-го порядка

$$\left\{ \begin{array}{l} y' = y_1\\ y'_1 = y_2\\ y'_2 = y_3\\ \ldots\\ y'_{(m-1)} = f(x,y_1,y_2,\ldots,y_{(m-1)} \end{array} \right.$$

Численным решением системы (12.9), а значит и дифференциального уравнения m-го порядка (12.8) является m табличных функций

$$y, y_1=y', y_2=y'', \ldots, y_m=y^{(m-1)},$$

т.е. функция y(x) и все ее производные, включая производную (m-1)-го порядка.

При этом каждая из табличных функций определяется на промежутке [a, b] с шагом h и включает n узловых точек. Таким образом, численным решением уравнения (12.8) или системы (12.9) является матрица порядка $$(n\times m)$$, (табл. 12.8)

i	x	y	y1=y',	y1=y''1,	:	y_m-1=y^(m-1)
0	x0	y0	(y1)0	(y2)0	:	(ym-1)0
1	x1	y1	(y1)1	(y2)1	:	(ym-1)1
2	x2	y2	(y1)2	(y2)2	:	(ym-1)2
3	x3	y3	(y1)3	(y2)3	:	(ym-1)3
:	:	:	:	:	:	:
n	xn	yn	(y1)n	(y2)n	:	(ym-1)n

где

m - порядок дифференциального уравнения, равен количеству столбцов матрицы,

n = (b-a)/h - количество шагов интегрирования, равно количеству строк матрицы.

Каждый j -й столбец матрицы - это массив решений одной j -й табличной функции по всем n шагам интегрирования.

Каждая i -ая строка матрицы - это массив решений m табличных функций на одном i -ом шаге интегрирования.

На графике решением дифференциального уравнения m -го порядка (12.8) является совокупность $$(n\times m)$$ узловых точек. При этом каждому шагу интегрирования, т.е. каждому значению xi, $$i=\overline{1,n}$$, соответствуют m узловых точек скоординатами

$$(x_i,y_1), (x_i,(y_1)_i), (x_i,(y_2)_i), \ldots, (x_i,(y_{m-1})_i).$$

При построении алгоритма задачи будем как и ранее расчет вести по шагам интегрирования, т.е. в цикле по $$i=\overline{1,n}.$$ При этом, как и ранее, на каждом i -ом шаге цикла будем рассчитывать решение дифференциального уравнения и тут же его печатать. Тогда нет необходимости формировать матрицу решений, а можно ограничиться формированием массива решений (12.15), который соответствует одной i-й строке матрицы:

$$Y=[y(1), y(2), y(3), \ldots , y(m)],$$

где

$$y(1)= y(x),\\ y(2)= y'(x),\\ y(3)= y"(x),\\ \ldots \ldots\\ y(m)= y^{m-1}(x).$$

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

Рассмотрим пример (12.19).

Дано дифференциальное уравнение второго порядка

$$y''=\frac{1}{x}y' - \frac{x^2-1}{x^2}y$$

С учетом обозначений (12.18) имеем:

$$y(1)=y(x),\\ y(2)=y'(x)$$

Тогда дифференциальное уравнение (12.19) сводится к системе:

$$\left\{ \begin{array}{l} y'(1)=y(2),\\ y'(2)=\frac{1}{x}y(2)-\frac{x^2-1}{x^2}y(1). \end{array} \right.$$

Вычисление правых частей уравнений системы (12.20) будем выполнять в подпрограмме PRAV. При этом подпрограмма PRAV будет иметь вид

$$\left\{ \begin{array}{l} f(1)=y(2),\\ f(2)=\frac{1}{x}y(2)-\frac{x^2-1}{x^2}y(1). \end{array} \right.$$

Обращение к подпрограмме

PRAV(m, x, Y, F),

где

m -порядок системы,

x - значение x на i -м шаге интегрирования,

Y=[y(1), y(2), y(3), : , y(m)] - входной массив длиной m,

F=[f(1), f(2),:, f(m)] - результат, массив значений правых частей уравнений системы (12.10) длиной m.

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

основной программы, в которой организуем циклический процесс по всем шагам интегрирования, $$i=\overline{1,n}$$ ;

подпрограммы RGK, в которой на каждом i -м шаге реализуется метод Рунге Кутта (4-го порядка) для системы дифференциальных уравнений m -го порядка. Здесь на каждом шаге интегрирования вычисляется массив решений Y длиной m. Для этого внутри подпрограммы RGK организуем циклы по j, где $$j=\overline{1,m}$$ ;

подпрограммы PRAV, обращение к которой осуществляется из подпрограммы RGK для вычисления массива F длиной m - значений правых частей уравнений системы.

Схема алгоритма основной программы представлена на рис.12.17.

(рис 12.17) Схема алгоритма основной программы

Здесь

m -порядок системы,

h -шаг интегрирования,

n -количество шагов интегрирования,

x -начальное и далее - текущее значение x,

Y -массив длиной m, куда заносим начальные и далее - текущие значения решений системы на одном шаге интегрирования.

В подпрограмме RGK для вычисления элементов массива Y, используем те же формулы, что и для решения одного дифференциального уравнения 1-го порядка методом Рунге-Кутта (4-го порядка), но с учетом поправки на массивы.

Тогда

$$T(1,j)=h\cdot f_ j(x, Y),\\ y1(j)=y(j)+T(1,j)/2, j=\overline{1,m},\\ T(2,j)=h\cdot f_ j(x+h/2, Y1),\\ y1(j)=y(j)+T(2,j)/2, j=\overline{1,m},\\ T(3,j)=h\cdot f_ j(x+h/2, Y1),\\ y1(j)=y(j)+T(3,j)/2, j=\overline{1,m},\\ T(4,j)=h\cdot f_ j(x+h, Y1),\\ y(j)=y(j)+\frac{1}{6}(T(1,j)+2T(2,j)+2T(3,j)+T(4,j)).$$

Здесь

Y - массив решений длиной m,

Y1 - рабочий массив длиной m,

T - рабочая матрица порядка ( $$4\times m$$ ).

Для вычисления значений fj правых частей дифференциальных уравнений системы перед вычислением каждого элемента матрицы Т на каждом j -ом шаге необходимо выполнить обращение к подпрограмме PRAV.

Схема алгоритма подпрограммы RGK представлена на рис.12.18

(рис 12.18) Схема алгоритма подпрограммы RGK