Здравствуйте.
Проблема в следующем. Мне нужно получить решение ОДУ в виде сплайна. Очевидный путь - получить значения функции в наборе точек посредством
ode45(..), сшить эти точки сплайном
spline(..).
Этот способ мне не нравится по 2 причинам. Для построения сплайна необходимо знать значения функции и производной в точках. Эти значения при решении ДУ мы, вообще говоря, рассчитываем. То есть нет никакого смысла вызывать функцию
spline(..). Во-первых, это оверхед (чёрт с ним), во вторых, что более важно, функция
spline(..) не использует информацию о значении производной функции (а ведь она очень полезна), а "угадывает" значение производной сама. В связи с этим увеличивается ошибка полученного решения.
Для пояснения своей мысли приведу пример (написан на scilab'e, нет большой разницы с matlab'ом - функции аналогичны):
////////////////////////////////////
// solve DE x' = -x^2
function [dx] = rhs(t, x)
dx = - x ^ 2;
endfunction
////////////////////////////////////
// solving...
t = [0:0.1:1]';
func_values = ode(1, 0, t, rhs)';
////////////////////////////////////
// take a spline
deriv_values = rhs(t, func_values);
deriv_values2 = splin(t, func_values);
function y = f1(q)
y = (interp(q, t, func_values, deriv_values) - 1 ./ (1 + q))^2;
endfunction
function y = f2(q)
y = (interp(q, t, func_values, deriv_values2) - 1 ./ (1 + q))^2;
endfunction
// calculate integral quadratic error
disp(intg(0, 1, f1)); // err = 1.648D-12
disp(intg(0, 1, f2)); // err = 9.577D-11
То есть если я вычисляю коэффициенты сплайна "вручную", то квадратичная ошибка получается на 2 порядка меньше, нежели при вызове функции
spline(..).
Собственно вопрос: как в данной ситуации лучше поступить? Умеет ли функция
ode45(..) возвращать помимо набора значений функции ещё и набор значений производной? Или же сразу структуру, описывающую сплайн? Или, может, помимо
ode45 существуют более удобные функции?
