Solve ODE using backward euler's method

x' = λ - ρx - βxz;
y' = βxz - δy;
z' = py - cz;
x0=43100; y0 = 0, z0 = 0.0033, λ = 388, ρ = 0.009 δ = 0.18, p = 50000, c = 23, β=3.61e-8

Is there a built-in function in matlab to solve the above non-linear system using the backward euler's method?

Answer 1

[Torsten]

Initialize

x_old = 43100, y_old = 0 and z_old = 0.0033

Compute x_new by solving the nonlinear system of equations

(x_new-x_old)/dt = lambda - rho*x_new - beta*x_new*z_new
(y_new-y_old)/dt = beta*x_new*z_new - delta*y_new
(z_new-z_old)/dt = p*y_new - c*z_new

by fixed-point iteration or with MATLAB's fsolve, e.g.

This gives you the solution for your system at time t=dt.

Set

x_old = x_new, y_old = y_new and z_old = z_new

and solve the above system again for x_new, y_new and z_new.

This gives you the solution at time t=2*dt.

Continue until you reach t=tfinal.

Best wishes

Torsten.

[Walter Roberson]

Torsten, that looks more like Forward Euler ??

[Torsten]

As far as I know, Forward Euler evaluates the right-hand side with the old values:

(x_new-x_old)/dt = lambda - rho*x_old - beta*x_old*z_old
(y_new-y_old)/dt = beta*x_old*z_old - delta*y_old
(z_new-z_old)/dt = p*y_old - c*z_old

Best wishes

Torsten.

[KC]

Thank you Mr. Torsten, I will try this code!

Answer 2

[Walter Roberson]

No, there is no built-in function in MATLAB for that.

[Walter Roberson]

But perhaps one of the File Exchange contributions would be useful: http://www.mathworks.com/matlabcentral/fileexchange/31427-1st-oder-ode-solver-euler-backward-method-based

[KC]

Is there an example somewhere of how to solve a system of ODE's using the backward euler's method? I would want to understand the concept first, so I can implement it in MATLAB. I googled for quite some time but was not able to find a proper example. Thanks!