CODEWITHSUNDEEP
mathmatlabnumericaldifferential-equations
SimpleWater
SimpleWater

Reputation: 113

How to obtain the numerical solution of these differential equations with matlab

I have differential equations derived from epidemic spreading. I want to obtain the numerical solutions. Here's the equations,

enter image description here

t is a independent variable and ranges from [0,100]. The initial value is 

y1 = 0.99; y2 = 0.01; y3 = 0;

At first, I planned to deal these with ode45 function in matlab, however, I don't know how to express the series and the combination. So I'm asking for help here.

**

The problem is how to express the right side of the equations as the odefun, which is a parameter in the ode45 function.

**

Upvotes: 1

Views: 775

Answers (2)

peter karasev
peter karasev

Reputation: 2614

 function demo(a_in)
 X = [0;0;0];
  T = [0:.1:100];
  a = a_in; % for nested scope

  [Xout,  Tout ]= ode45( @myFunc, T, X );

    function [dxdt] = myFunc( t, x )
         % nested function accesses "a"
         dxdt = 0*x + a; 
         % Todo: real value of dxdt. 
    end
  end

What about this, and you simply need to fill in the dxdt from your math above? It remains to be seen if the numerical roundoff matters...

Edit: there's a serious issue due to the 1=y1+y2+y3 constraint. Is that even allowed, since you have an IVP with 3 initial values given and 3 first order ODE's? If that constraint is a natural consequence of the equations, it may not be needed.

Upvotes: 1

Michael J. Barber
Michael J. Barber

Reputation: 25052

Matlab has functions to calculate binomial coefficients (number of combinations) and the finite series can be expressed just as matrix multiplication. I'll demonstrate how that works for the sum in the first equation. Note the use of the element-wise "dotted" forms of the arithmetic operators.

Calculate a row vector coefs with the constant coefficients in the sum as:

octave-3.0.0:33> a = 0:20;
octave-3.0.0:34> coefs = log2(a * 0.05 + 1) .* bincoeff(20, a);

The variables get combined into another vector:

octave-3.0.0:35> y1 = 0.99;
octave-3.0.0:36> y2 = 0.01;
octave-3.0.0:37> z = (y2 .^ a) .* ((1 - y2) .^ a) .* (y1 .^ a);

And the sum is then just evaluated as the inner product:

octave-3.0.0:38> coefs * z'

The other sums are similar.

Upvotes: 2

Related Questions