MATLAB: Solve ODE using backward euler’s method

backward euler's

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?

Best Answer

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.

Related Solutions

MATLAB: 4th order Runge-Kutta code that can solve for several intial conditions

Here's a version:

% Initial conditions
t0=0;
x0=0:0.2:1;
y0=1:0.2:2;
dt=0.1;
tz=1;
t_range= t0:dt:tz;
X = zeros(numel(x0),numel(t_range));
Y = zeros(numel(y0),numel(t_range));
x_rk=zeros(1,numel(t_range));
y_rk= zeros(1,numel(t_range));
for n = 1:numel(x0)
    x_rk(1)=x0(n);
    y_rk(1)=y0(n);
    for i=1:length(t_range)-1
        [x_rk(i+1), y_rk(i+1)]=rk_method(x_rk(i),y_rk(i),dt);
    end
    X(n,:) = x_rk;
    Y(n,:) = y_rk;
    %% PLots
    subplot(numel(x0)/2,2,n)
    plot (t_range, x_rk,'b-o',t_range, y_rk,'y-o','MarkerSize',5);
    grid
    title([num2str(x0(n)),'   ',num2str(y0(n))])
    legend('x','y')
end
disp('X = '), disp(X)
disp('Y = '), disp(Y)
%% Functions
function [dxdt, dydt]= f(x,y)
dxdt =  5*x-3*y;
dydt = -6*x+2*y;
end 
function [x1, y1]=rk_method(x0,y0,dt)
[k11, k12]=f(x0,y0);     
[k21, k22]=f(x0+0.5*dt*k11,y0+0.5*dt*k12); 
[k31, k32]=f(x0+0.5*dt*k21,y0+0.5*dt*k22); 
[k41, k42]=f(x0+dt*k31, y0+dt*k32); 
x1=x0+dt*((k11/6)+((k21+k31)/3)+(k41/6));
y1=y0+dt*((k12/6)+((k22+k32)/3)+(k42/6));
end

MATLAB: Changing a parameter in a differential equation

Hi Sean,

you can pass the value of epsilon into the ode function with

[tSol,YSol]=ode45(@(t,Y) quantumsystem(t,Y,epsilon),tRange,Y0);

and

function dYdt = quantumsystem(t,Y,epsilon);
    x=Y(1);
    y=Y(2);
    z=Y(3);
    gamma=1/2;
    c=1;
    dtheta=gamma*c/(t^2 + c^2);
    dxdt=(-1/epsilon)*y - dtheta*z; 
    dydt=(1/epsilon)*x;
    dzdt=dtheta*x;
    dYdt=[dxdt;dydt;dzdt];
end

Since ode45 picks the time increments, different values of epsilon may lead to different lengths of tSol and YSol, but the information is there. One of the easier ways to deal with that issue is to set up a for loop with different values of epsilion and use

plot(tSol,z)
hold on

inside the loop as you are doing, only without the color specification.

Best Answer

Related Solutions

MATLAB: 4th order Runge-Kutta code that can solve for several intial conditions

MATLAB: Changing a parameter in a differential equation

Related Question