Coupled systems of ODE's solved by matlabs ODE45 (Runge Kutta Alg.)

In summary, this function takes in various parameters and computes the right hand sides of the state and sensitivity equations for an age-structured omnivory model. It outputs a column vector of dy/dt values, including the adult and juvenile predator and consumer densities, resource density, and sensitivities for each species. The parameters used in the model include eRP, eCP, eRC, aRP, aCP, aRC, hRP, hCP, hRC, mP, mC, nP, r, K, and mJ. These parameters are used to calculate the partial derivatives and update the state and sensitivity values.
  • #1
DrWahoo
53
0
Use the code as needed, its free to all!
Please note this code is only for one mfile, you will need to add them all to the same file path for it to work. Others will be in the next post. They have to run simultaneously to work and give output.
Code: 
function dydt = rhs_both_age(~,y,flag,p,eRP,eCP,eRC,aRP,aCP,aRC,hRP,hCP,hRC,mP,mC,nP,r,K,mJ)
%--------------------------------------------------------------------------
% Simultaniously computes the right hand sides of the state and sensitivity 
% equations for the age structure omnivory model.
%
% Input:    t = time. Not used, but needed by ode45
%           flag = (not used) placeholder for compatibility with ode45 y =
%           column vector (length 6) of states and sensitivities
%               y(1) = adult predator density, y(2) = juvenile predator
%               density y(3) = consumer density, y(4) = resource density, 
%               y(5) = adult predator sensitivity, y(6) = juvenile predator
%               sensitivity y(7) = consumer sensitivity y(8) = resource
%               sensitivity
%           p = integer that tells which parameter sensitivity to compute
%           eRP, eCP, ... , r, K are model parameters
%
% Output:   dydt = column vector (length 8) of dy/dt values
%--------------------------------------------------------------------------

%---the column vector of four zeros on top and the partials of f with 
%   respect to the model parameters. 
if p == 1 %eRP
    partial = [0;
                0;
                0;
                0;
                0;
                (aRP*y(4)*y(1))/(1+aRP*hRP*y(4)+aCP*hCP*y(3));
                0;
                0]; 
elseif p == 2 %eCP
    partial = [0;
                0;
                0;
                0;
                0;
                (aCP*y(3)*y(1))/(1+aRP*hRP*y(4)+aCP*hCP*y(3));
                0;
                0]; 
elseif p == 3 %eRC
    partial = [0;
                0;
                0;
                0;
                0;
                0;
                (aRC*y(4)*y(3))/(1+aRC*hRC*y(4));
                0]; 
elseif p == 4 %aRP
    partial = [0;
                0;
                0;
                0;
                0;
                (eRP*y(4)*y(1)*(1+aCP*hCP*y(3))-eCP*aCP*hRP*y(4)*y(3)*y(1))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2);
                (aCP*hRP*y(4)*y(3)*y(1))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2);
                -(y(4)*y(2))/((1+aRP*hRP*y(4))^2)-(y(4)*y(1)*(1+aCP*hCP*y(3)))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2)];  
elseif p == 5 %aCP
    partial = [0;
                0;
                0;
                0;
                0;
                (eCP*y(3)*y(1)*(1+aRP*hRP*y(4))-eRP*aRP*hCP*y(4)*y(3)*y(1))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2);
                -(y(3)*y(1)*(1+aRP*hRP*y(4)))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2);
                (aRP*hCP*y(4)*y(3)*y(1))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2)]; 
elseif p == 6 %aRC
    partial = [0;
                0;
                0;
                0;
                0;
                0;
                (eRC*y(4)*y(3))/((1+aRC*hRC*y(4))^2);
                -(y(4)*y(3))/((1+aRC*hRC*y(4))^2)]; 
elseif p == 7 %hRP
    partial = [0;
                0;
                0;
                0;
                0;
                -(aRP*y(4)*y(1)*(eRP*aRP*y(4)+eCP*aCP*y(3)))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2);
                (aRP*aCP*y(4)*y(3)*y(1))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2);
                (((aRP)^2)*((y(4))^2)*y(2))/((1+aRP*hRP*y(4))^2)+(((aRP)^2)*((y(4))^2)*y(1))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2)];  
elseif p == 8 %hCP
    partial = [0;
                0;
                0;
                0;
                0;
                -(aCP*y(3)*y(1)*(eRP*aRP*y(4)+eCP*aCP*y(3)))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2);
                (((aCP)^2)*((y(3))^2)*y(1))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2);
                (aRP*aCP*y(4)*y(3)*y(1))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2)]; 
elseif p == 9 %hRC
    partial = [0;
                0;
                0;
                0;
                0;
                0;
                -(eRC*((aRC)^2)*((y(4))^2)*y(3))/((1+aRC*hRC*y(4))^2);
                (((aRC)^2)*((y(4))^2)*y(3))/((1+aRC*hRC*y(4))^2)];             
                
elseif p == 10 %mPA
    partial = [0;
                0;
                0;
                0;
                -y(1);
                0;
                0;
                0]; 
elseif p == 11 %mC
    partial = [0;
                0;
                0;
                0;
                0;
                0;
                -y(3);
                0]; 
elseif p == 12 %nP
    partial = [0;
                0;
                0;
                0;
                y(2);
                -y(2);
                0;
                0];             
elseif p == 13 %r
    partial = [0;
                0;
                0;
                0;
                0;
                0;
                0;
                y(4)*(1-(y(4)/K))];  
elseif p == 14 %K
    partial = [0;
                0;
                0;
                0;
                0;
                0;
                0;
                ((r*(y(4))^2)/(K)^2)];   
 elseif p == 15 %mPJ
     partial = [0;
                0;
                0;
                0;
                0;
                -y(2);
                0;
                0];            
            
end

dydt_temp = [
nP*y(2)-mP*y(1);

((eRP*aRP*y(4)+eCP*aCP*y(3))/(1+aRP*hRP*y(4)+aCP*hCP*y(3)))*y(1)-(nP+mJ)*y(2);

y(3)*((eRC*aRC*y(4))/(1+aRC*hRC*y(4))-(aCP*y(1))/(1+aRP*hRP*y(4)+aCP*hCP*y(3))-mC);

y(4)*(r*(1-(y(4)/K))-(aRC*y(3))/(1+aRC*hRC*y(4))-(aRP*y(2))/(1+aRP*hRP*y(4))-(aRP*y(1))/(1+aRP*hRP*y(4)+aCP*hCP*y(3)));

-mP*y(5)+nP*y(6);

((eRP*aRP*y(4)+eCP*aCP*y(3))/(1+aRP*hRP*y(4)+aCP*hCP*y(3)))*y(5)-(nP+mP)*y(6)+((eCP*aCP*y(1)*(1+aRP*hRP*y(4))-eRP*aRP*aCP*hCP*y(4)*y(1))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2))*y(7)+((eRP*aRP*y(1)*(1+aCP*hCP*y(3))-eCP*aRP*aCP*hRP*y(3)*y(1))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2))*y(8);

-((aCP*y(3))/(1+aRP*hRP*y(4)+aCP*hCP*y(3)))*y(5)+((eRC*aRC*y(4))/(1+aRC*hRC*y(4))-(aCP*y(1)*(1+aRP*hRP*y(4)))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2)-mC)*y(7)+((eRC*aRC*y(3))/((1+aRC*hRC*y(4))^2)+(aCP*aRP*hRP*y(3)*y(1))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2))*y(8);

-((aRP*y(4))/(1+aRP*hRP*y(4)+aCP*hCP*y(3)))*y(5)-((aRP*y(4))/(1+aRP*hRP*y(4)))*y(6)-((aRC*y(4))/(1+aRC*hRC*y(4))+(aRP*aCP*hCP*y(4)*y(1))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2))*y(7)+(r*(1-((2*y(4))/K))-(aRC*y(3))/((1+aRC*hRC*y(4))^2)-(aRP*y(2))/((1+aRP*hRP*y(4))^2)-(aRP*y(1)*(1+aCP*hCP*y(3)))/((1+aRP*hRP*y(4)+aCP*hCP*y(3))^2))*y(8)

];

dydt = dydt_temp + partial;
 

FAQ: Coupled systems of ODE's solved by matlabs ODE45 (Runge Kutta Alg.)

1. What is a coupled system of ODEs?

A coupled system of ODEs is a set of differential equations that are connected or dependent on each other. This means that the solution of one equation affects the solution of the other equations in the system.

2. How does MATLAB's ODE45 solve coupled systems of ODEs?

MATLAB's ODE45 uses the Runge-Kutta algorithm to solve coupled systems of ODEs. This algorithm is a numerical method that uses a combination of function evaluations at different points to approximate the solution of the differential equations.

3. What is the advantage of using MATLAB's ODE45 to solve coupled systems of ODEs?

MATLAB's ODE45 is a highly efficient and accurate method for solving coupled systems of ODEs. It can handle a wide range of differential equations and provides flexible options for controlling the accuracy of the solution. Additionally, it has built-in error control mechanisms to ensure reliable results.

4. How do I specify the initial conditions and parameters for a coupled system of ODEs in MATLAB's ODE45?

To solve a coupled system of ODEs using ODE45, you need to first define the differential equations in a function file. Then, you can specify the initial conditions and parameters in the function call when using ODE45. Alternatively, you can also use the "odefun" input argument to specify the function directly in the ODE45 function call.

5. Can I use ODE45 to solve a system of ODEs with a varying number of equations?

Yes, ODE45 can handle systems of ODEs with a varying number of equations. You can define the differential equations in a function file and pass the appropriate number of initial conditions and parameters in the function call. Alternatively, you can use the "odefun" input argument to specify the function directly in the ODE45 function call.

Similar threads

MHB ODE 45 solving coupled ODE's (Pred, Prey, Resource) with graph outputs.
Replies
1
Views
2K
MATLAB How can I plot a Hohmann Transfer Orbit in MATLAB using ode45?
Replies
1
Views
2K
MATLAB A Lorenz's system of ODEs doesn't get executed in Matlab
Replies
5
Views
2K
MATLAB Animating Vector Addition in Matlab (Looking for improvements)
Replies
4
Views
1K
MATLAB Question about how solve 2 body problem in orbit using MATLAB ode45
Replies
6
Views
3K
MATLAB Problem with Generating Barnsley Fern Fractal in MATLAB
Replies
3
Views
3K
Using ODE45 on Matlab to Integrate A Problem
Replies
1
Views
1K
MATLAB How Do You Solve Coupled ODEs Using MATLAB's ode45?
Replies
1
Views
2K
MATLAB Plotting the tragectory of an asteroid in MATLAB
Replies
7
Views
3K
MATLAB Projectile motion with friction in MATLAB (ODE45)
Replies
4
Views
7K

Hot Threads

Recent Insights

Back
Top