Advanced Search

Question: SI Model Simulation: Error in Jacobian Matrix Calculation and Eigenvalue Computation

Question:SI Model Simulation: Error in Jacobian Matrix Calculation and Eigenvalue Computation

Paras31 90 Maple 2024
matrix syntax linearalgebra diff eigenvalues
May 02 2024
1

SI_MODEL.mw
I have tried the following code to simulate the following SI Model. But when I run it, it gives me several mistakes and most importantly does not calculate the Jacobian Matrix and the eigenvalues. 

Any help would be deeply appreciated


Many of the most interesting dynamics in nature have to do with interactions between organisms. These interactions are often subtle, indirect, and difficult to detect. Interactions in which one organism consumes all or part of another. This includes predator-prey, herbivore-plant, and parasite-host interactions. These linkages are the prime movers of energy through food chains. They are an important factor in the ecology of populations, determining the mortality of prey and the birth of new predators.
​
Mathematical models and logic suggest that a coupled system of predator and prey should cycle: predators increase when prey is abundant, prey are driven to low numbers by predation, the predators decline, and the prey recover, ad infinitum. One such model that simulates predator-prey interactions is the Lotka - Volterra Model.

 


We will discuss a behavior of 2 D - subsystems: SI - model.


SI Model (without predator)

 

We will study a mathematical model that appears in eco - eco-epidemiology.


The SI model describes the interactions between S-prey with density x and I-prey with density y under the assumption z≡0


Break the population into compartments.

• 

Prey (S)
• 

Infected prey (I)


SI Equations :

(D(x))(t) = r*x(t)-b*x(t)^2-c*x(t)*y(t)-beta*x(t)*y[t]/(a+x(t))

(D(y))(t) = -mu*y(t)+beta*x(t)*y(t)/(a+x(t))

– 

 r stands for the intrinsic growth rate of S-prey.
– 

b defines the intra-class competition in S-prey.
– 

c characterizes the inter-class competition between S-prey and I-prey.
– 

μ stands for the mortality of the I-prey
– 

β considered as the bifurcation parameter.
• 

Note: The dynamics of the system depend on parameter β

 

Modeling with Differential Equations

 

• 

We will fix the parameters in the study to have the following values r = 1, b = 1, c = 1/100, μ = 4/10, a=1/2.75.
• 

The initial values will be x(0)=0.2, y(0)=0.05


As mentioned above the dynamics of the system depend on parameter β. So we will consider two cases when β=0.75 and β=1


We will solve the system for the above values of parameters and initial conditions and for a time interval [0,200] , then we will create its plots

 

Case 1

  restart

beta := .75; r := 1; b := 1; c := 1/100; mu := 4*(1/10); a := 1/2.75; f := r*x(t)-b*x(t)^2-c*x(t)*y(t)-(3/4)*x(t)*y(t)/(a+x(t)); g := -mu*y(t)+(3/4)*x(t)*y(t)/(a+x(t)); deq1 := diff(x(t), t) = f; deq2 := diff(y(t), t) = g; equilibrio := solve({f = 0, g = 0}, {x(t), y(t)}); jacobian := simplify(Matrix(2, 2, [[diff(f, x(t)), diff(f, y(t))], [diff(g, x(t)), diff(g, y(t))]])); A := simplify(eval(jacobian, equilibrio[1])); B := simplify(eval(jacobian, equilibrio[2])); Q := simplify(eval(jacobian, equilibrio[3])); eigen1 := Eigenvalues(A); eigen2 := Eigenvalues(B); eigen3 := Eigenvalues(Q); soln := dsolve({deq1, deq2, x(0) = .2, y(0) = 0.5e-1}, numeric, output = listprocedure); plots:-display(plots:-odeplot(soln, [t, x(t)], 0 .. 200, color = blue, legend = ["x(t)"]), plots:-odeplot(soln, [t, y(t)], 0 .. 200, color = red, legend = ["y(t)"]), labels = ["t", "Population"], title = "Population Dynamics")

  

NULL

 

 

Case 2

restart``

beta := 1; r := 1; b := 1; c := 1/100; Mu := 4*(1/10); a := 1/2.75; f := r*x(t)-b*x(t)^2-c*x(t)*y(t)-x(t)*y(t)/(a+x(t)); g := -Mu*y(t)+x(t)*y(t)/(a+x(t)); deq1 := diff(x(t), t) = f; deq2 := diff(y(t), t) = g; equilibrio := solve({f = 0, g = 0}, {x(t), y(t)}); jacobian := Matrix(2, 2, [[diff(f, x(t)), diff(f, y(t))], [diff(g, x(t)), diff(g, y(t))]], simplify); J_at_equilibrio := [eval(jacobian, equilibrio)]; eigen1 := Eigenvalues(A); eigen2 := Eigenvalues(B); eigen3 := Eigenvalues(Q); soln := dsolve({deq1, deq2, x(0) = .2, y(0) = 0.5e-1}, numeric); plots:-display(plots:-odeplot(soln, [t, x(t)], 0 .. 200, color = blue, legend = ["x(t)"]), plots:-odeplot(soln, [t, y(t)], 0 .. 200, color = red, legend = ["y(t)"]), labels = ["t", "Population"], title = "Population Dynamics with Beta = 1")

  

NULL


 

Download SI_MODEL.mw

 

 
Please Wait...