MATLAB хочет преобразовать явный алгоритм Эйлера в неявный алгоритм Эйлера для СИСТЕМЫ ODE 1-го порядка - PullRequest
Новая
       50

MATLAB хочет преобразовать явный алгоритм Эйлера в неявный алгоритм Эйлера для СИСТЕМЫ ODE 1-го порядка

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

Прежде всего я не могу использовать ЛЮБЫЕ встроенные решатели ODE для этого.Я закодировал эту систему ODE явным методом Эйлера, но вместо этого мне нужно переписать ее неявным методом Эйлера.Если я просто переключаю «i» на «i + 1» (например, «y (1, i)» на «y (1, i + 1)»), то ответ получается неприлично неправильным.

A_init= 4; 
B_init= 1.1;
C_init= 4;
y0 = [A_init; B_init; C_init];
h = 20/100;
n = 100;
t0 = 0;
tf = 20;
t = linspace(t0,tf,n);
y = zeros(numel(y0) , n); 
y(:, 1) = y0(:);

dAdt= @(a,b) 7.27*(b-a*b+ a-(8.375*10^-5)*a^3);
dBdt= @(b,a,c)(-b-a*b+c)/77.27;
dCdt= @(c,a) 0.4*(a-c);

for i = 1:n-1
    y(1,i+1) = y(1,i) +h*feval(dAdt,y(1,i),y(2,i));
    y(3,i+1) = y(3,i) +h*feval(dCdt,y(3,i),y(1,i));
    y(2,i+1) = y(2,i) +h*feval(dBdt,y(2,i),y(1,i),y(3,i));
end

1 Ответ

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

Ваша идея верна, но помните, что если вы измените i в правой части на i+1, у вас будет y(1,i+1), y(2,i+1), y(3,i+1), появляющийся в обеих сторонах уравнения.Это означает, что вы на самом деле должны решить для y(1,i+1), y(2,i+1) и y(3,i+1) на каждом шаге.Поскольку ваши три уравнения связаны, вы должны решать нелинейную систему уравнений на каждом временном шаге, используя fsolve или fzero.

Читатьответ на этот вопрос , который показывает, как это сделать для случая с одним уравнением.

...