0

Before anything, I can't use ANY built-in ODE solvers for this. I coded this system of ODEs with explicit euler's method but I need to instead rewrite it with implicit euler's. If I just switch the "i" to "i+1" (ex "y(1,i)" to "y(1,i+1)") then the answer comes out obscenely wrong. 

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
Dcg
  • 1

1 Answers1

0

Your idea is correct, but remember that if you change the i on the right-hand side to i+1, you will have y(1,i+1), y(2,i+1), y(3,i+1) appearing on both sides of the equation. This means that you actually have to solve for y(1,i+1), y(2,i+1) and y(3,i+1) at each step. Since your three equations are coupled, you have to solve a nonlinear system of equations at each time-step, using fsolve or fzero.

Read the answer to this question, which shows how to do this for the single equation case.

Savithru
  • 735
  • 6
  • 12