Программное обеспечение. Ниже приведен текст подпрограммы на алгоритмическом языке FORTRAN

Ниже приведен текст подпрограммы на алгоритмическом языке 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,...) ,

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

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

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

Поиск по сайту: