ode

Аргументы

y0: вектор или матрица вещественных значений, начальные условия.
t0: вещественный скаляр, начальное время.
t: вещественный вектор, моменты времени для которых находится решение.
f: внешняя функция (функция, строка или список), правая сторона дифференциального уравнения.
type: строка, тип используемой программы решения. Имеются следующие типы программ решения: "adams", "stiff", "rk", "rkf", "fix", "discrete", "root".
atol, rtol: absolute and relative tolerances on the final solution y (decimal numbers). If each is a single value, it applies to each component of y. Otherwise, it must be a vector of same size as size(y), and is applied element-wise to y.
jac: внешняя функция (функция, строка или список), якобиан функции f.
ng: целое число.
g: внешняя функция (функция, символьная строка или список).
k0: целое число (начальное время).
kvect: целочисленный вектор.
y: вещественный вектор или матрица (ВЫХОДНАЯ).
rd: вещественный вектор (ВЫХОДНОЙ).
w, iw: вещественные векторы (ВХОДНЫЕ/ВЫХОДНЫЕ). See ode() optional output

Описание

ode решает явные обыкновенные дифференциальные уравнения, определённые как:

$\begin{eqnarray} \begin{array}{l} \dfrac{dy}{dt} = f(t,y)\\ y(t_0)=y_0 \end{array} \end{eqnarray}$

Это интерфейс для различных программ решения, в частности для ODEPACK.

В данной справке описывается использование ode для стандартных явных систем ОДУ.

Самый простой способ вызова ode: y = ode(y0, t0, t, f), где y0 - вектор начальных условий, t0 - начальное время, t - вектор моментов времени, для которых вычисляется решение y и y - матрица векторов решения y=[y(t(1)),y(t(2)),...].

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

Если f является Scilab-функцией, то её последовательность вызова должна быть:
```
ydot = f(t,y)
```
где t - вещественный скаляр (время), а y - вещественный вектор (состояние) и ydot - вещественный вектор (производная первого порядка dy/dt).
Если f - строка, то это - имя компилированной процедуры Fortran или функции C. Например, если вызывается ode(y0, t0, t, "fex"), то вызывается подпрограмма fex.

Процедура Fortran должна иметь заголовок:
```
fex(n,t,y,ydot)
```
где n - целое число, t -- скаляр двойной точности, y и ydot - векторы двойной точности.

Функция C должна иметь заголовок:
```
fex(int *n,double *t,double *y,double *ydot)
```
где t - время, y - состояние, а ydot - производная состояния (dy/dt).

Эта внешняя функция может быть собрана способом, независящим от ОС, используя ilib_for_link и динамически связана с Scilab с помощью функции link.
Может случиться, что симулятору f потребуются дополнительные аргументы. В этом случае можно использовать следующую возможность. Аргумент f может также быть списком lst=list(simuf,u1,u2,...un), где simuf является Scilab-функцией с синтаксисом: ydot = f(t,y,u1,u2,...,un), а u1, u2, ..., un - дополнительные аргументы, автоматически передаваемые симулятору simuf.

Функция f может возвращать вместо вектора матрицуp x q. С этой матричной нотацией решается система n=p+q ОДУ dY/dt=F(t,Y), где Y - матрица p x q. Затем начальные условия, Y0, должны быть так же матрицей p x q, а результат ode - матрицей p x q(T+1) [Y(t_0),Y(t_1),...,Y(t_T)].

Допуски rtol и atol являются порогами для относительной и абсолютной оценённых ошибок. Оценённая ошибка y(i) равна: rtol(i)*abs(y(i))+atol(i) и интегрирование выполняется пока эта ошибка мала для каждого из элементов состояния. Если rtol и/или atol является константой, то rtol(i) и/или atol(i) являются набором констант. Значения по умолчанию для rtol и atol соответственно равны rtol=1.d-5 и atol=1.d-7 для большинства программ решения, а для "rfk" и "fix" rtol=1.d-3 и atol=1.d-4.

Для жёстких задач лучше указывать якобиан функции правой стороны в качестве необязательного аргумента jac. Якобиан является внешней функцией, т. е. функцией со специальным синтаксисом, либо именем процедуры Fortran или функции C (символьная строка) с определённой последовательностью вызова, либо списком.

Если jac является функцией, то синтаксис должен быть J=jac(t,y), где t -- вещественный скаляр (время), а y - вещественный вектор (состояние). Результирующая матрица J должна вычислять df/dx, т. е. J(k,i) = dfk/dxi, где fk - k-тый элемент f.
Если jac является символьной строкой, то она ссылается на имя процедуры Fortran или функции C.

Процедура Fortran должна иметь заголовок:
```
subroutine fex(n,t,y,ml,mu,J,nrpd)
integer n,ml,mu,nrpd
double precision t,y(*),J(*)
```
Функция C должна иметь заголовок:
```
void fex(int *n,double *t,double *y,int *ml,int *mu,double *J,int *nrpd,)
```
В большинстве случаев не нужно ссылаться на ml, mu и nrpd.
Если jac является списком, то для f применяются те же договорённости.

Необязательные аргументы w и iw являются векторами для хранения информации, возвращаемой подпрограммой интегрирования (см. ode_optional_output). Когда эти векторы указываются с правой стороны ode, интегрирование начинается заново с теми же параметрами, что и до остановки.

Программам решения ODEPACK можно передать больше опций с помощью переменной %ODEOPTIONS. See odeoptions.

Программы решения

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

<не задано>:: Программа решения lsoda из пакета ODEPACK вызывается по умолчанию. Она автоматически выбирает между нежёстким методом прогноза-исправления Адамса (predictor-corrector Adams method) и жёстким методом обратной дифференциальной формулой (ОДФ) ( Backward Differentiation Formula (BDF) method). Изначально она применяет нежёсткий метод и динамически проверяет данные для того, чтобы решить какой метод использовать.
"adams":: Используется для нежёстких задач. Вызывается программа решения lsode из пакета ODEPACK, и она использует метод Адамса.
"stiff":: Это для жёстких задач. Вызывается программа решения lsode из пакета ODEPACK, и она использует метод ОДФ.
"rk":: Адаптивный метод Рунге-Кутты 4-го порядка (RK4).
"rkf":: Используется программа Шампайна и Уотса (Shampine и Watts), основанная на методе Рунге-Кутты-Фелберга пары 4 и 5-го порядка (RKF45). Она для нежёстких и среднежёстких задач, когда получаемые вычисления не дороги. Этот метод не следует, в общем случае, использовать, если пользователь требует высокую точность.
"fix":: Та же программа решения, что и "rkf", но пользовательский интерфейс очень прост, т. е. программе решения могут быть переданы только параметры rtol и atol. Это самый простой метод для пробы.
"root":: Программа решения ОДУ с возможностью нахождения корней. Используется программа решения lsodar из пакета ODEPACK. Она является вариантом программы решения lsoda, где она ищет корни заданной векторной функции. См. справку по ode_root.
"discrete":: Моделирование дискретного времени. См. справку по ode_discrete.

Примеры

В следующем примере мы решим обыкновенное дифференциальное уравнение dy/dt=y^2-y*sin(t)+cos(t) с начальным условием y(0)=0. Используем программу решения по умолчанию.

function ydot=f(t, y)
    ydot = y^2 - y*sin(t) + cos(t)
endfunction
y0 = 0;
t0 = 0;
t = 0:0.1:%pi;
y = ode(y0,t0,t,f);
plot(t,y)

В следующем примере мы решаем уравнение dy/dt=A*y. Точное решение равно y(t)=expm(A*t)*y(0), где expm - матричная экспонента. Неизвестное равно матрице y(t) размером 2 на 1.

function ydot=f(t, y)
    ydot = A*y
endfunction
function J=Jacobian(t, y)
    J = A
endfunction
A = [10,0;0,-1];
y0 = [0;1];
t0 = 0;
t = 1;
ode("stiff",y0,t0,t,f,Jacobian)
// Сравним с точным решением:
expm(A*t)*y0

В следующем примере мы решаем ОДУ dx/dt = A*x(t) + B*u(t) с u(t)=sin(omega*t). Обратите внимание, что дополнительные аргументы f: A, u, B, omega передаются в f в списке.

function xdot=linear(t, x, A, u, B, omega)
    xdot = A*x + B*u(t,omega)
endfunction
function ut=u(t, omega)
    ut = sin(omega*t)
endfunction
A = [1 1;0 2];
B = [1;1];
omega = 5;
y0 = [1;0];
t0 = 0;
t = [0.1,0.2,0.5,1];
ode(y0,t0,t,list(linear,A,u,B,omega))

В следующем примере мы решим дифференциальное уравнение Риккати dX/dt=A'*X + X*A - X'*B*X + C, где начальное условие X(0) является единичной матрицей 2 на 2.

function Xdot=ric(t, X, A, B, C)
    Xdot = A'*X + X*A - X'*B*X + C
endfunction
A = [1,1;0,2];
B = [1,0;0,1];
C = [1,0;0,1];
y0 = eye(A);
t0 = 0;
t = 0:0.1:%pi;
X = ode(y0,t0,t,list(ric,A,B,C))

В следующем примере мы решаме дифференциальное уравнение dY/dt=A*Y, где неизвестная Y(t) является матрицей 2 на 2. Точное решение равно expm Y(t)=expm(A*t), где expm -- матричная экспонента.

function ydot=f(t, y, A)
    ydot = A*y;
endfunction
A = [1,1;0,2];
y0 = eye(A);
t0 = 0;
t = 1;
ode(y0,t0,t,list(f,A))
// Сравнение с точным решением
expm(A*t)
ode("adams",y0,t0,t,f)

С компилятором

Следующий пример требует компилятор C.

// ---------- Простое одномерное ОДУ (внешняя функция на языке C)
ccode=['#include <math.h>'
       'void myode(int *n,double *t,double *y,double *ydot)'
       '{'
       '  ydot[0]=y[0]*y[0]-y[0]*sin(*t)+cos(*t);'
       '}']
mputl(ccode,TMPDIR+'/myode.c') // создаём C-файл
// Компилируем
cd TMPDIR
ilib_for_link('myode','myode.c',[],'c',[],'loader.sce');
exec('loader.sce') // пошаговая компоновка
y0 = 0;
t0 = 0;
t = 0:0.1:%pi;
y = ode(y0,t0,t,'myode');

Литература

Alan C. Hindmarsh, "lsode and lsodi, two new initial value ordinary differential equation solvers", ACM-Signum newsletter, vol. 15, no. 4 (1980), pp. 10-11.

Используемые функции

Связанные процедуры могут быть найдены в директории SCI/modules/differential_equations/src/fortran: lsode.f lsoda.f lsodar.f

Смотрите также

odeoptions — установка опций для программ решения ОДУ
ode_optional_output — описание необязательных выходных данных программы решения ОДУ
ode_root — программа решения обыкновенных дифференциальных уравнений с поиском корней
ode_discrete — программа решения обыкновенных дифференциальных уравнений, моделирование дискретного времени
dae — Differential algebraic equations solver
impl — дифференциальное алгебраическое уравнение
odedc — программа решения дискретно-непрерывных ОДУ
csim — simulation (time response) of linear system
ltitr — discrete time response (state space)
rtitr — discrete time response (transfer matrix)
intg — определённый интеграл

Report an issue
<< numderivative	Дифференциальные уравнения	ode_discrete >>

Синтаксис