Библиотека функций MATLAB’а

 

Важным средством, позволяющим ускорить разработку программ, являются библиотеки стандартных программ. Есть такая библиотека и в MATLAB’е.

Надо научиться находить в HELP’е описания нужных программ. В частности в конце этого занятия вам будет предложено написать программу решения дифференциальных уравнений по методу Эйлера. Написав свою программу, вам надо будет найти функцию MATLAB’а, решающую эту же задачу и используя ее сравнить решения полученные вашей программой и библиотечной

 

Что еще почитать

 

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

Надо сказать, что такой литературы, которая учила бы принципам программирования, а не синтаксическим правилам какого-либо языка программирования не так уж много. «Программистов часто учат программировать и редко – отладке» (ван Тассел).

Могу порекомендовать следующие книги, в которых вопросы отладки и тестирования подробно рассмотрены:

1. Ван Тассел Д. Стиль, разработка, эффективность, отладка и испытания программ. 1981
2. Керниган, Пайк. Практика программирования. 2004
3. Иванова Г.С. Технология программирования. 2002

 

 

Задания

 

Написать программу решения дифференциальных уравнений по методу Эйлера. При отладке и тестировании вести учет ошибок. Получить решение этой же задачи, используя библиотечную функцию MATLAB’а. Сравнить полученные результаты

 

№	Задача Коши для уравнения первого порядка
Задача Коши для уравнения высшего порядка

 

 

Справочная информация

Обыкновенное дифференциальное уравнение 1-го порядка

связывает независимую переменную t,искомую функцию y и ее производную. Решение дифференциальногоуравнениязаключается в отысканиифункции y = y(t), обращающей это уравнение в тождество на конечном или бесконечном интервале (a, b). Различают общее и частные решения дифференциального уравнения. Общее решение имеет вид y = y(t, C), где C - произвольная постоянная интегрирования. Его графическим отображением является семейство кривых (см. рис.1), называемых интегральными. Каждая интегральная кривая является отображением частного решения, соответствующего своему значению постоянной C. Для выделения частного решения из множества общего решения необходимо задать начальное условие

.

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

Для решения задачи Коши существует множество методов, которые делятся на одношаговые и многошаговые. Все они позволяют получить искомое решение в виде таблично заданной функции, в той или иной мере согласующееся с истинным частным решением (см. рис.2). Эти группы методов различаются объемом информации, которая используется для вычисления координат очередной точки табличной функции. Одношаговые методы используют значения функции и ее производной только в одной предыдущей точке, в то время как многошаговые - в нескольких. К одношаговым методам решения задачи Коши относятся метод Эйлера, модифицированный метод Эйлера, методы Рунге-Кутта и другие.

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

,

или в матричной форме следующим образом

, ,

где

, , .

Система дифференциальных уравнений связывает независимую переменную x, искомые функции y ₁, y ₂ ,..., y_n и их первые производные. В данном случае решение задачи Коши заключается в отыскании функций y ₁ = y ₁ (x), y ₂ = y ₂ (x),..., y_n = y_n(x), обращающих каждое уравнение системы в тождество на конечном или бесконечном интервале (a, b) и удовлетворяющих начальным условиям.

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

,

, , , …, ,

то замена

, , , …, ,

приводит ее к

,

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

Для решения такой задачи Коши используются те же методы, что для обыкновенных дифференциальных уравнений 1-го порядка. Это обуславливается тем, что матричная форма записи задачи Коши для нормальной системы полностью совпадает с ее формулировкой для этих уравнений. Единственным отличием является то, что вместо функций y(t) и f(t,y) используются вектор-функции y и f, состоящие из n функций y ₁ (t), y ₂ (t),..., y_n(t) и f ₁ (t,y ₁,..., y_n), f ₂ (t,y ₁,..., y_n),..., f_n(t,y ₁,..., y_n), соответственно. При этом расчетные схемы методов и оценки их погрешностей сохраняются.

 

 

Метод Эйлера

 

Он является старейшим методом решения задачи Коши и заключается в последовательном применении следующих формул

,

,

,

геометрическая интерпретация которых представлена на рис.3. В точке t ₀ вычисляется значение производной dy/dt через f(t,y), которое определяет угол наклона касательной к графику точного решения задачи Коши. Следующая точка численного решения определяется как точка на этой касательной с абсциссой t ₁ = t ₀ + h. В компактном виде эти соотношения записываются следующим образом

,

.

Метод Эйлера - самый простой метод численного интегрирования. Он относится к методам первого порядка точности, поскольку его решение совпадает с истинным только в том случае, когда истинное решение является линейной функцией y = a ₁ + a ₂ t. Его погрешность e_k на каждом шаге пропорциональна величине h ². Это обусловлено тем, что в качестве направления, определяющего положение следующей точки численного решения, используется касательная в крайней левой точке каждого отрезка [t _k ,t _k+1 ]. Из рис.3 видно, что для получения более точного численного решения надо использовать некоторое промежуточное направление между направлениями касательных в крайних точках рассматриваемого отрезка.

Соотношения метода Эйлера для нормальной системы в матричной форме имеют вид

,

.

или в развернутой форме

, .

Геометрическая интерпретация работы метода Эйлера решения задачи Коши для нормальной системы идентична его геометрической интерпретации для дифференциальных уравнений 1-го порядка. Однако в данном случае движение осуществляется вдоль некоторойгиперкривойв (n+1) -мерном пространстве переменных t, y ₁, y ₂ ,...,y_n.

 

 

Программное обеспечение

 

Ниже приведен текст подпрограммы на алгоритмическом языке FORTRAN, реализующей метод Эйлера:

Subroutine Eu(n, dt, k, Fun, t, Y)

 

В качестве параметров в подпрограмме используются:

n - порядок системы дифференциальных уравнений;

dt - длина отрезка интегрирования уравнений;

k - количество шагов на длине отрезка интегрирования;

Fun - имя внешней подпрограммы типа Subroutine, с помощью которой вычисляются значения вектора правой части f (t, y)системы;

t - значение аргумента системы в начале отрезка интегрирования при обращении к подпрограммам и в конце этого отрезка (t + dt)после отработки подпрограмм;

Y - массив (n элементов) значений вектора решения y (t)системы уравнений в начале отрезка интегрирования при обращении к подпрограммам и в конце этого отрезка после отработки подпрограмм.

Подпрограмма Fun оформляется в виде

 

Subroutine Fun(n, t, Y, F)

Dimension Y(n), F(n)

F(1) =...........

..................

F(n) =...........

End

 

Здесь F - массив значений вектора правой части f (t, y) системы уравнений.

 

Subroutine Eu(n, dt, k, Fun, t, Y)

Real F1(n)

h = dt/k

Do m = 1, k

Call Fun(n, t, Y, F1)

t = t + h

Do i = 1, n

Y(i) = Y(i) + h*F1(i)

End do

End do

 

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

· [t,Y] = ode23(fun, tspan, y0) – решает задачу Коши для системы обыкновенных дифференциальных уравнений с использованием комбинации методов Рунге-Кутта относительно невысокого порядка (2-го и 3-го);

· [t,Y] = ode45(fun, tspan, y0) - решает задачу Коши для системы обыкновенных дифференциальных уравнений с использованием комбинации методов Рунге-Кутта более высокого порядка (4-го и 5-го).

Входные аргументы указанных функций:

fun - функция для вычисления вектора правой части системы, интерфейс которой должен иметь вид:

F = fun(t, Y),

где t – независимый параметр системы, Y- вектор значений неизвестных функций, а выходной параметр F – вектор вычисленных компонент правой части;

tspan - вектор, определяющий интервал интегрирования системы. Если этот вектор состоит из двух компонент, то их значения задают начало и конец интервала. Если требуется, чтобы в качестве решения были выданы значения искомых функций в конкретно заданных точках интервала, то число компонент должно быть не менее трех (можно задавать в виде имени массива или при помощи конструктора массивов, например в виде [2:0.1:5]);

y0 - вектор начальных значений.

Выходные аргументы:

t - вектор, состоящий из значений независимой переменной, которым соответствуют вычисленные значения решения, помещенные в массив Y (если вектор tspan состоит лишь из двух компонент, то количество значений автоматически определяется функцией с целью наиболее точного отображения на графике);

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

Параметр fun, являющийся фактическим параметром, указывающим предварительно написанную функцию для вычисления правой части системы, задается в виде:

@fname, где fname - имя m-функции.

Точность получаемого решения по умолчанию регулируется условием обеспечения относительной погрешности не более 10^-3 или абсолютной погрешности не более 10^-6. Если требуется задать иные требования точности, то функции следует вызывать с дополнительным аргументом options:

[t,Y] = ode23(fun, tspan, y0, options) или

[t,Y] = ode45(fun, tspan, y0, options).

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

options = odeset('name1',value1,'name2',value2,...),

где name i - название управляющего параметра, а value i - новое присваиваемое ему значение. Среди управляющих параметров имеются RelTol и AbsTol, задающие относительную и абсолютную погрешности решения систем. Так, например, если предполагается, что вызываемая функция должна обеспечить относительную погрешность решения не более 10^-5 или абсолютную не более 10^-8, то следует выполнить оператор

options = odeset('RelTol',1e-5,'AbsTol',1e-8),

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