Calculate ODE's with MATLAB

Question

The differential equations:

α'(t)=s(β-βα+α-qα^2)

β'(t)=(s^-1)(-β-αβ+γ)

γ'(t)=w(α-γ)

Intitial values

α (0) = 30.00

β (0) = 1.000

γ (0) = 30.00

Calculation

I want to solve the problem from t_0=0 to t=10 while using the values s = 1, q = 1 and w = 0.1610

I've no idea how to write the function for the ODE's and would very much appreciate the help!

Answer 1

I'm not usually in the habit of solving other people's homework, but today's your lucky day I guess.

So, you have a system of coupled ordinary differential equations:

α'(t) = s(β-α(β+1)-qα²)

β'(t) = (-β-αβ+γ)/s

γ'(t) = w(α-γ)

and you want to solve for

y = [α(t) β(t) γ(t)]

with 0 < t < 10, s = 1, q = 1, w = 0.1610. The way to do that in matlab is to define a function that computes the derivative ([α'(t) β'(t) γ'(t)]), and throw this in one of the ODE solvers (ode45 is a good first bet):

s = 1;
q = 1; 
w = 0.1610;

% Define y(t) = [α(t) β(t) γ(t)] = [y(1) y(2) y(3)]:

deriv = @(t,y) [...
    s * (y(2) - y(1)*(y(2)+1) - q*y(1)^2)   % α'(t)
    (-y(2) - y(1)*y(2) + y(3))/s            % β'(t)
    w * (y(1)-y(3))                         % γ'(t)
];

% initial value
y0 = [30 1 10];

% time span to integrate over
tspan = [0 10];

% solve ODE numerically
[t, y] = ode45(deriv, tspan, y0)

This will output

y =
   30.0000    1.0000   10.0000
   28.5635    0.9689   10.0049   % numerical solutions at times t
   27.2558    0.9413   10.0094
   26.0603    0.9166   10.0136
   ...        ...      ...
   = α(t)     = β(t)   = γ(t)

t =
         0
    0.0016
    0.0031   % corresponding times
    0.0047
    ...

You can plot all this like so:

figure, clf, hold on
plot(t, y(:,1), 'r')
plot(t, y(:,2), 'g')
plot(t, y(:,3), 'b')
legend('\alpha(t)', '\beta(t)', '\gamma(t)')

which results in this figure: