Неправильное значение для производной решения ОДУ в MATLAB - PullRequest
0 голосов
/ 11 декабря 2018

Попытка решить задачу Коши для нелинейного ОДУ с помощью классического метода Рунге-Кутты, описанного в здесь (см. ode4).После выполнения этой части

q = 10;
T = 5;
N = @(g1, g2) g1^2 + sin(g2);

t = linspace(0, T, q);
GCC = nonlinearGreenCC(N, 1, T);
di = diff(GCC);

я оцениваю GCC(1) и di(1).В то время как GCC(1) = 0 как ожидалось, di(1) = 1.6e-05.Не могу понять почему, потому что условие Коши на первой производной равно 1. Как исправить неточность / ошибку?Любая помощь высоко ценится.

Функция nonlinearGreenCC выглядит следующим образом:

function G = nonlinearGreenCC(N, a0, T)

h = .0001;
eqGreen = @(t, g)[g(2); - N(g(1), g(2))];
Cc = [0, a0];
G = ode4(eqGreen, 0, h, T, Cc);

end

1 Ответ

0 голосов
/ 11 декабря 2018

GCC - это массив пар значений, начинающийся для более точной интеграции с чем-то вроде

[  0.00000000e+00,   1.00000000e+00],
[  9.99957927e-05,   9.99915855e-01],
[  1.99983171e-04,   9.99831715e-01],
[  2.99962136e-04,   9.99747579e-01],
[  3.99932687e-04,   9.99663448e-01],

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

 0                         1
 1.66652642851402e-05      1.00001666526429
-1.40231141152723e-05      0.999985976885885
 1.66638620438405e-05      1.00001666386204
-1.40217119017434e-05      0.999985978288098

Не содержит выборки временного интервала, действительно, на уровне, который вы называете diff, значение h неизвестно.diff не может вычислить коэффициенты разницы, которые вы ожидаете.Если вы прочитаете документацию, вы обнаружите, что она вычисляет только различия.И первое различие возвращает только первое значение 1.66652642851402e-05.


Прямая ошибка, которая заставляет алгоритм ode4 давать странные результаты, состоит в том, что eqGreen создает векторы столбцов, где он должен возвращать векторы строк,Поскольку начальное значение представляет собой вектор строки, результат добавления векторов строки и столбца в ode4 создает матрицу 2x2, которая добавляется к результату в виде 2 строк с непонятными последствиями.Наличие обоих векторов строк помещает результаты в один ряд с переменным значением и производной.Измените на

eqGreen = @(t, g)[g(2), - N(g(1), g(2))];
Cc = [0, a0];

, и результат будет

GCC(1:5,:) =

                     0                     1
  9.99957927208439e-05     0.999915855174488
  0.000199983171186392     0.999831714893813
  0.000299962135851046     0.999747579156327
  0.000399932687169042     0.999663447960381

di(1:5,:)=

  9.99957927208439e-05  -8.41448255122224e-05
  9.99873784655483e-05  -8.41402806746050e-05
  9.99789646646540e-05  -8.41357374861129e-05
  9.99705513179961e-05  -8.41311959463020e-05
  9.99621384254097e-05  -8.41266560547282e-05

Если вы масштабируете последнее путем деления на h=1e-4, вы получите ожидаемые результаты.

...