Решение ОДУ 2 степени в Matlab модифицированным методом Эйле

[Imaginarium]

Здравствуйте.

Прошу помощи в разборе задачи: нужно решить ОДУ 2 степени в матлаб, вида
$y'' = F(y,t)$, только вручную, модифицированным (с пересчетом) методом Эйлера.

Для начала нужно сделать обыкновенный гармонический осциллятор, чтобы убедиться, что все понимаю правильно. Самое главное- вынести правую часть в отдельный m-файл. И то ли лыжи не едут, то ли как обычно все, но гармонический осциллятор не получается.
Уравнение использую вида $y'' = -y$
Правую часть пытаюсь сделать сразу, на будущее, зависящей от $y$ и $t$.
Вот код:
файл euler_s.m

Используется синтаксис Matlab M

clc;
clear all;
 
sec = 1;
h = 0.2; %шаг интегрирования
 
zp = zeros(n, sec/h+1); %predictor 
zc = zeros(n, sec/h+1); %corrector 2st func.
 
yp = zeros(n, sec/h+1); %predictor 
yc = ones(n, sec/h+1); %corrector 1st func.
 
for t = 2:sec/h+1
    %prediction
    yp(t) = yc(t-1) + h*zc(t-1);
    zp(t) = zc(t-1) + h*Fun(zc(t-1),t-1);
    %correction
    yc(t) = yc(t-1)+h*(zc(t-1) + zp(t))/2 
    zc(t) = zc(t-1)+h*(Fun(zc(t-1),t-1) + Fun(zp(t),t))/2;
    
end
 
plot([1:sec/h+1],yc),grid;
 

Здесь уравнение 2 степени с помощью известной замены становится системой ОДУ 1 степени:
$ \begin{cases} y' = z\\ z' = Fun(y,t) \end{cases} $

Далее файл правой части Fun.m:

Используется синтаксис Matlab M

function = Fun(r,t);
Fun = -r;
 

Вот правая часть как раз для гармонического осциллятора. График показывает горизонтальную прямую на высоте начального условия. Не понимаю ошибки у себя, прошу помощи.

Спасибо.

[ewert]

Imaginarium в сообщении #846634 писал(а):

zp = zeros(n, sec/h+1); %predictor
zc = zeros(n, sec/h+1); %corrector 2st func.

У Вас вроде бы n нигде не инициализировано. И в любом случае организация вычислений какая-то нелепая. Во-первых: зачем массив-то двумерный?... Во-вторых, нет решительно никакой необходимости заводить массив для хранения результатов прогноза -- он нужен только для результатов коррекции.

-- Вт апр 08, 2014 12:31:57 --

Да, и ещё: раз уж на будущее -- то следует предусмотреть зависимость правой части не только от времени и от игрека, но и от производной этого игрека.