Return to Answer

corrected code

edited Mar 27, 2022 at 15:00

2.7k
6
19

 (*   set up constants and ODEs  *)  Subscript[\[Mu], B] = 0.5;   Subscript[\[Beta], 1] = 0.1;  Subscript[\[Beta], 3] = 0.2;  Subscript[\[Gamma], 1] = 0.7;  Subscript[\[Gamma], 2] = 0.221;  Subscript[\[Phi], b1] = 0.9;  Subscript[\[Phi], w1] = 0.05;  Subscript[\[Phi], w2] = 0.08;  Subscript[\[Beta], 2] = 0.8;  Subscript[\[Beta], 4] = 0.2;  Subscript[\[Mu], C] = 0.5;   Subscript[\[Mu], E] = 0.01;  Subscript[\[Phi], b2] = 0.6;  Subscript[\[Gamma], 3] = 0.109;  Subscript[\[Mu], D] = .07;  Subscript[\[Gamma], 4] = 0.219;  ode1 = Derivative[1][b][t] == -(b[t] Subscript[\[Mu], B]) +   b[t] Subscript[\[Beta],   1] (-b[t] - c[t] + 1) + (1 - Subscript[\[Beta],   1]) Subscript[\[Beta], 3] d[t] (-b[t] - c[t] + 1) -   b[t] Subscript[\[Gamma], 1] - b[t] e[t] Subscript[\[Phi], b1] -   b[t] c[t] Subscript[\[Phi], w1];  ode2 = Derivative[1][c][t] == Subscript[\[Beta], 2]  Subscript[\[Beta], c[t]2]c[t] (-b[t] - c[t] + 1) + (1 - Subscript[\[Beta], 2]) Subscript[\[Beta], 4] e[t] (-b[t] - c[t] + 1) +   b[t] e[t] Subscript[\[Phi], b1] +   b[t] c[t] Subscript[\[Phi], w1] - Subscript[\[Gamma], 2] c[t] -   c[t] Subscript[\[Mu], C];  ode3 = Derivative[1][d][t] ==   b[t] Subscript[\[Beta],   1] (1 - Subscript[\[Beta], 3]) (-d[t] - e[t] + 1) -   d[t] c[t] Subscript[\[Phi], b2] - Subscript[\[Gamma], 3] d[t] -   d[t] Subscript[\[Mu], D] +   Subscript[\[Beta], 3] d[t] (-d[t] - e[t] + 1) -   d[t] e[t] Subscript[\[Phi], w2];  ode4 = Derivative[1][e][t] ==   d[t] c[t] Subscript[\[Phi], b2] +   Subscript[\[Beta],   2] (1 - Subscript[\[Beta], 4]) c[t] (-d[t] - e[t] + 1) +   Subscript[\[Beta], 4] e[t] (-d[t] - e[t] + 1) +   d[t] e[t] Subscript[\[Phi], w2] - Subscript[\[Gamma], 4] e[t] -   e[t] Subscript[\[Mu], E];  (*   run system with some arbitrary initial conditions  *)  {fb, fc, fd, fe} =   NDSolveValue[{ode1, ode2, ode3, ode4, b[0] == 0.5, c[0] == -0.2,   d[0] == 1.1, e[0] == 0.5}, {b, c, d, e}, {t, 0, 1000}];  (*   plot each function separately  *)    Plot[{fb[t], fc[t], fd[t], fe[t]}, {t, 0, 1000},   PlotLabel -> Style["fb[t],fc[t],fd[t],fe[t]", 16]]  (*   choose {fb[t],fc[t],fd[t]} and plot a 3D phase portrait of these \  fuctions  *)  ParametricPlot3D[{fb[t], fc[t], fd[t]}, {t, 0, 100},   BoxRatios -> {1, 1, 1},   PlotLabel -> Style["Phase portrait of {fb[t],fc[t],fd[t]}", 16]]

(* set up constants and ODEs *) Subscript[\[Mu], B] = 0.5; Subscript[\[Beta], 1] = 0.1; Subscript[\[Beta], 3] = 0.2; Subscript[\[Gamma], 1] = 0.7; Subscript[\[Gamma], 2] = 0.221; Subscript[\[Phi], b1] = 0.9; Subscript[\[Phi], w1] = 0.05; Subscript[\[Phi], w2] = 0.08; Subscript[\[Beta], 2] = 0.8; Subscript[\[Beta], 4] = 0.2; Subscript[\[Mu], C] = 0.5; Subscript[\[Mu], E] = 0.01; Subscript[\[Phi], b2] = 0.6; Subscript[\[Gamma], 3] = 0.109; Subscript[\[Mu], D] = .07; Subscript[\[Gamma], 4] = 0.219; ode1 = Derivative[1][b][t] == -(b[t] Subscript[\[Mu], B]) + b[t] Subscript[\[Beta], 1] (-b[t] - c[t] + 1) + (1 - Subscript[\[Beta], 1]) Subscript[\[Beta], 3] d[t] (-b[t] - c[t] + 1) - b[t] Subscript[\[Gamma], 1] - b[t] e[t] Subscript[\[Phi], b1] - b[t] c[t] Subscript[\[Phi], w1]; ode2 = Derivative[1][c][t] == Subscript[\[Beta], 2]  c[t] (-b[t] - c[t] + 1) + (1 - Subscript[\[Beta], 2]) Subscript[\[Beta], 4] e[t] (-b[t] - c[t] + 1) + b[t] e[t] Subscript[\[Phi], b1] + b[t] c[t] Subscript[\[Phi], w1] - Subscript[\[Gamma], 2] c[t] - c[t] Subscript[\[Mu], C]; ode3 = Derivative[1][d][t] == b[t] Subscript[\[Beta], 1] (1 - Subscript[\[Beta], 3]) (-d[t] - e[t] + 1) - d[t] c[t] Subscript[\[Phi], b2] - Subscript[\[Gamma], 3] d[t] - d[t] Subscript[\[Mu], D] + Subscript[\[Beta], 3] d[t] (-d[t] - e[t] + 1) - d[t] e[t] Subscript[\[Phi], w2]; ode4 = Derivative[1][e][t] == d[t] c[t] Subscript[\[Phi], b2] + Subscript[\[Beta], 2] (1 - Subscript[\[Beta], 4]) c[t] (-d[t] - e[t] + 1) + Subscript[\[Beta], 4] e[t] (-d[t] - e[t] + 1) + d[t] e[t] Subscript[\[Phi], w2] - Subscript[\[Gamma], 4] e[t] - e[t] Subscript[\[Mu], E]; (* run system with some arbitrary initial conditions *) {fb, fc, fd, fe} = NDSolveValue[{ode1, ode2, ode3, ode4, b[0] == 0.5, c[0] == -0.2, d[0] == 1.1, e[0] == 0.5}, {b, c, d, e}, {t, 0, 1000}]; (* plot each function separately *) Plot[{fb[t], fc[t], fd[t], fe[t]}, {t, 0, 1000}, PlotLabel -> Style["fb[t],fc[t],fd[t],fe[t]", 16]] (* choose {fb[t],fc[t],fd[t]} and plot a 3D phase portrait of these \ fuctions *) ParametricPlot3D[{fb[t], fc[t], fd[t]}, {t, 0, 100}, BoxRatios -> {1, 1, 1}, PlotLabel -> Style["Phase portrait of {fb[t],fc[t],fd[t]}", 16]]

 (*   set up constants and ODEs  *)  Subscript[\[Mu], B] = 0.5;   Subscript[\[Beta], 1] = 0.1;  Subscript[\[Beta], 3] = 0.2;  Subscript[\[Gamma], 1] = 0.7;  Subscript[\[Gamma], 2] = 0.221;  Subscript[\[Phi], b1] = 0.9;  Subscript[\[Phi], w1] = 0.05;  Subscript[\[Phi], w2] = 0.08;  Subscript[\[Beta], 2] = 0.8;  Subscript[\[Beta], 4] = 0.2;  Subscript[\[Mu], C] = 0.5;   Subscript[\[Mu], E] = 0.01;  Subscript[\[Phi], b2] = 0.6;  Subscript[\[Gamma], 3] = 0.109;  Subscript[\[Mu], D] = .07;  Subscript[\[Gamma], 4] = 0.219;  ode1 = Derivative[1][b][t] == -(b[t] Subscript[\[Mu], B]) +   b[t] Subscript[\[Beta],   1] (-b[t] - c[t] + 1) + (1 - Subscript[\[Beta],   1]) Subscript[\[Beta], 3] d[t] (-b[t] - c[t] + 1) -   b[t] Subscript[\[Gamma], 1] - b[t] e[t] Subscript[\[Phi], b1] -   b[t] c[t] Subscript[\[Phi], w1];  ode2 = Derivative[1][c][t] == Subscript[\[Beta], 2]c[t] (-b[t] - c[t] + 1) + (1 - Subscript[\[Beta],2]) Subscript[\[Beta], 4] e[t] (-b[t] - c[t] + 1) +   b[t] e[t] Subscript[\[Phi], b1] +   b[t] c[t] Subscript[\[Phi], w1] - Subscript[\[Gamma], 2] c[t] -   c[t] Subscript[\[Mu], C];  ode3 = Derivative[1][d][t] ==   b[t] Subscript[\[Beta],   1] (1 - Subscript[\[Beta], 3]) (-d[t] - e[t] + 1) -   d[t] c[t] Subscript[\[Phi], b2] - Subscript[\[Gamma], 3] d[t] -   d[t] Subscript[\[Mu], D] +   Subscript[\[Beta], 3] d[t] (-d[t] - e[t] + 1) -   d[t] e[t] Subscript[\[Phi], w2];  ode4 = Derivative[1][e][t] ==   d[t] c[t] Subscript[\[Phi], b2] +   Subscript[\[Beta],   2] (1 - Subscript[\[Beta], 4]) c[t] (-d[t] - e[t] + 1) +   Subscript[\[Beta], 4] e[t] (-d[t] - e[t] + 1) +   d[t] e[t] Subscript[\[Phi], w2] - Subscript[\[Gamma], 4] e[t] -   e[t] Subscript[\[Mu], E];  (*   run system with some arbitrary initial conditions  *)  {fb, fc, fd, fe} =   NDSolveValue[{ode1, ode2, ode3, ode4, b[0] == 0.5, c[0] == -0.2,   d[0] == 1.1, e[0] == 0.5}, {b, c, d, e}, {t, 0, 1000}];  (*   plot each function separately  *)    Plot[{fb[t], fc[t], fd[t], fe[t]}, {t, 0, 1000},   PlotLabel -> Style["fb[t],fc[t],fd[t],fe[t]", 16]]  (*   choose {fb[t],fc[t],fd[t]} and plot a 3D phase portrait of these \  fuctions  *)  ParametricPlot3D[{fb[t], fc[t], fd[t]}, {t, 0, 100},   BoxRatios -> {1, 1, 1},   PlotLabel -> Style["Phase portrait of {fb[t],fc[t],fd[t]}", 16]]

Source Link

answered Mar 27, 2022 at 14:36

josh

2.7k
6
19

It's a good try pyring. But in order to solve first-order ODEs numerically, you need two things: (1) Supply numeric values to all constants, (2) Supply numeric values for initial conditions. You have 16 constants and four equations so that's 20 values you need to assign. In the code below I chose random values between 0 and 1. The important thing is to just get it running without syntax errors or run-time errors. Later you can adjust it to your needs. In the code below I plotted individually the four functions written as fb[t], fc[t],fd[t] and fe[t]. Then I chose three of them fb[t], fc[t] and fd[t] to construct a 3D phase portrait which appears as a simple line. Maybe with other choices of initial conditions and constants you can achieve something more interesting. Just cut this code and paste it into your notebook and it should run without errors. Then you can experiment with changing the initial conditions and constants.

(* set up constants and ODEs *) Subscript[\[Mu], B] = 0.5; Subscript[\[Beta], 1] = 0.1; Subscript[\[Beta], 3] = 0.2; Subscript[\[Gamma], 1] = 0.7; Subscript[\[Gamma], 2] = 0.221; Subscript[\[Phi], b1] = 0.9; Subscript[\[Phi], w1] = 0.05; Subscript[\[Phi], w2] = 0.08; Subscript[\[Beta], 2] = 0.8; Subscript[\[Beta], 4] = 0.2; Subscript[\[Mu], C] = 0.5; Subscript[\[Mu], E] = 0.01; Subscript[\[Phi], b2] = 0.6; Subscript[\[Gamma], 3] = 0.109; Subscript[\[Mu], D] = .07; Subscript[\[Gamma], 4] = 0.219; ode1 = Derivative[1][b][t] == -(b[t] Subscript[\[Mu], B]) + b[t] Subscript[\[Beta], 1] (-b[t] - c[t] + 1) + (1 - Subscript[\[Beta], 1]) Subscript[\[Beta], 3] d[t] (-b[t] - c[t] + 1) - b[t] Subscript[\[Gamma], 1] - b[t] e[t] Subscript[\[Phi], b1] - b[t] c[t] Subscript[\[Phi], w1]; ode2 = Derivative[1][c][t] == Subscript[\[Beta], 2] c[t] (-b[t] - c[t] + 1) + (1 - Subscript[\[Beta], 2]) Subscript[\[Beta], 4] e[t] (-b[t] - c[t] + 1) + b[t] e[t] Subscript[\[Phi], b1] + b[t] c[t] Subscript[\[Phi], w1] - Subscript[\[Gamma], 2] c[t] - c[t] Subscript[\[Mu], C]; ode3 = Derivative[1][d][t] == b[t] Subscript[\[Beta], 1] (1 - Subscript[\[Beta], 3]) (-d[t] - e[t] + 1) - d[t] c[t] Subscript[\[Phi], b2] - Subscript[\[Gamma], 3] d[t] - d[t] Subscript[\[Mu], D] + Subscript[\[Beta], 3] d[t] (-d[t] - e[t] + 1) - d[t] e[t] Subscript[\[Phi], w2]; ode4 = Derivative[1][e][t] == d[t] c[t] Subscript[\[Phi], b2] + Subscript[\[Beta], 2] (1 - Subscript[\[Beta], 4]) c[t] (-d[t] - e[t] + 1) + Subscript[\[Beta], 4] e[t] (-d[t] - e[t] + 1) + d[t] e[t] Subscript[\[Phi], w2] - Subscript[\[Gamma], 4] e[t] - e[t] Subscript[\[Mu], E]; (* run system with some arbitrary initial conditions *) {fb, fc, fd, fe} = NDSolveValue[{ode1, ode2, ode3, ode4, b[0] == 0.5, c[0] == -0.2, d[0] == 1.1, e[0] == 0.5}, {b, c, d, e}, {t, 0, 1000}]; (* plot each function separately *) Plot[{fb[t], fc[t], fd[t], fe[t]}, {t, 0, 1000}, PlotLabel -> Style["fb[t],fc[t],fd[t],fe[t]", 16]] (* choose {fb[t],fc[t],fd[t]} and plot a 3D phase portrait of these \ fuctions *) ParametricPlot3D[{fb[t], fc[t], fd[t]}, {t, 0, 100}, BoxRatios -> {1, 1, 1}, PlotLabel -> Style["Phase portrait of {fb[t],fc[t],fd[t]}", 16]]