MATLAB: Solve a system of equations

differential equationsequationsmodelode45solve

So I have written a system of equations and used ode45 to solve it. I was just wondering if there is a more efficient way to do it. I am creating an ODE model and will later use certain methods to find the unknown parameters, but for now I am just guessing random values.

I have written my model as

 function dydt = ODE_Model(t,y,p_1,p_2,p_3,d_1,d_2,d_3,a,b,Tau)
dydt = zeros(3,1); 
dydt(1) = p_1*(a*cos(t*2*pi/365-Tau)+b)-d_1*y(1); 
dydt(2) = p_2*(a*cos(t*2*pi/365-Tau)+b)-d_2*y(2); 
dydt(3) = p_3*y(1)-d_3*y(2)*y(3); 
end

and then solved it using

    tspan = [0 5];
y0 = [4 5 6];
p_1 = 1;
p_2 = 1;
p_3 = 1;
d_1 = 1;
d_2 = 1;
d_3 = 1;
a=1;
b=1;
Tau=1;
SolExp = ode45(@(t,y) ODE_Model(t,y,p_1,p_2,p_3,d_1,d_2,d_3,a,b,Tau), tspan, y0);
TimePoints = 0:0.1:2;
y = deval(SolExp, TimePoints);
clf
plot(TimePoints, y)

I wanted to do P and D as a combined vector, so i would have P(1) .. P(6) and then i could just write something like P = [1 1 1 1 1 1], but reading about ode45 it seems I can only use one variable? I am not sure. Is there a neater way to do what I want to do, or is this as neat it gets?

Thanks

Best Answer

These may provide you with an approach appropriate to your problem:

Monod kinetics and curve fitting (link)
Solving Coupled Differential Equations (link)
Parameter Estimation for a System of Differential Equations (link)

All of these involve estimating the parameters of a system of differential equations, using lsqcurvefit to do the regressions.

Related Solutions

MATLAB: Error in ode45 (Index exceeds matrix dimensions)

First, ‘tu’ has to be an (Nx2) matrix.

Second, do not use global variables! They make your code much more difficult to troubleshoot. Pass ‘p’ and ‘tu’ as extra parameters instead.

Try this:

function dx=myeqn(t,x,p,tu)
output=interp1(tu(:,1),tu(:,2),t);
dx(1,1)= (p(1)/p(2))*(AF_function_in(t)-x(1))- (p(4)/p(2))*(x(1)-x(2));
dx(2,1)=  (p(4)/p(3)) *(x(1)-x(2));
end

and then call it in your ODE solver as:

[t,x] = ode45(@(t,x) myeqn(t,x,p,tu), tspan, x0);

with ‘x0’ being your (2x1) initial conditions vector.

MATLAB: Odearguments error for fitting data with system of ODEs

The ‘dydt’ definition should be:

dydt = [0; 0];

that would define a column vector. The differential equation runs correctrly with that change.

The problem is that you are fitting one column off data (‘f’), however your differential equation returns two columns. This throws:

Error using lsqcurvefit (line 269)
Function value and YDATA sizes are not equal.

That can be solved by using ‘t’ instead of ‘tspan’:

[T,y] = ode45(@MSolveFit,t,conds);

This works:

D = load('fittingdata_t_f.txt');
t = D(:,1);
f = D(:,2);
pguess = rand(4,1)*1E+6;
p = lsqcurvefit(@MSolve,pguess,t,f);
    function y = MSolve(p,t)
        cond1 = 0;
        cond2 = 0;
        conds = [cond1, cond2];
       
        tspan = [t(1) t(end)];
       
        % y(1) = P(t), y(2) = M(t)
        [T,y] = ode45(@MSolveFit,t,conds);
                
        y = y(:,2);
       
        function dydt = MSolveFit(t,y)
            nc = 2; % # monomers for primary nucleation
            m0 = 10; % initial concentration of tau
            mt = m0 - y(2);
           
            % parameters :
            % p(1) = kn
            % p(2) = km
            % p(3) = ka
            % p(4) = kp
           
            % dydt(1) = dPdt(t)
            % dydt(2) = dMdt(T)
            dydt = [0; 0];
            dydt(1) = p(1).*mt.^nc + p(2).*y(2) - p(3).*y(1).^2;
            dydt(2) = p(4).*mt.*y(1);          
        end
    end
figure
plot(t, f, 'p')
hold on
plot(t, MSolve(p,t), '-r')
hold off
grid

However I cannot get a good fit to your data, probably because I have no idea what you are doing or what the initial parammeter estimates should be.

Related Question