首先,我不能使用任何内置的ODE解算器来解决这个问题。我用显式的欧拉方法编写了这个系统的常微分方程,但是我需要用隐式的欧拉方法重写它。如果我只是把"i“改成"i+1”(例如"y(1,i)“to "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
https://stackoverflow.com/questions/53587066
复制