0
$\begingroup$

I have 42 coupled differential equations (6 of them are PDE's). My system has 6 objects and an object have 7 coupled equations. Since these objects are interacting with each other, their coupled equations also coupled with each other.

I can solve the system for 14 equations (2 objects). However, I can not run the code for 6 objects, NDSolve doesn't even give an error or warning. This is my first trial, and I am trying to make the simplest possible agent based model (i guess it can be referred as agent based modelling).

Is it possible to solve such a system with MMA? If so how can it be done? My final question, is there a different terminology for systems with high number of equations?

Any help would be appreciated.

The code:

Clear["Global`*"]
(*SET UP*)
tmax = 2500;
(*INITIAL CONDITIONS*)

inical = SetPrecision[ 0.200165185171545551`5. , 5];
iniope = 1/(1 + inical^2);
\[Mu]\[Mu] = SetPrecision[0.5500000, 5]
\[Mu]0 = SetPrecision[0.5155877999949406`, 5]
(*PARAMETERS*)
Dc = 20;
kflux = 8.1;
l = 20;(*--------*)
b = 0.111 // Rationalize;
k1 = 0.7;
\[Gamma] = 2.0;
k\[Gamma] = 0.1;
\[Beta] =(*0 *)0.02;
\[Tau]n = 2.0;
k2 = 0.7;
K1 = ((kflux \[Tau]n)/k1);
\[CapitalGamma]\[CapitalGamma] = ((\[Tau]n \[Gamma] 0.9)/k1);
\[CapitalGamma] = ((\[Tau]n \[Gamma] )/k1);
T = (\[Gamma]*\[Tau]n)/k1;
K = (k\[Gamma]/k1);
K2 = k2/k1;
B = \[Beta]/k1;
D0 = Dc \[Tau]n/l^2;
Dp0 = Dp \[Tau]n /l^2;
Dp = 300;
\[Mu]1 = 0.433;
 kp = 0.02;
(*q=-3.4;*)
xxx = -0.1;
qqq = -0.4;
Q1 = 0.09;
Q2 = 0.5;
ttt = 0.6;
zzz = -1.2;
sol = NDSolve[{
  D[v1o1[t, x], t] ==  
   D0 D[v1o1[t, x], x, 
      x] +((\[Mu]0) K1 ((b + v1o1[t, x])/(1 + v1o1[t, x])) v2o1[t, 
         x]) - ((\[CapitalGamma]\[CapitalGamma] v1o1[t, x])/(K + 
       v1o1[t, x])) + 
    B -10^zzz(*Q110^xxx/\[Tau]n*)(v7o1[t, x]) v1o1[t, x],
  D[v2o1[t, x], t] == 1/(1 + v1o1[t, x]^2) - v2o1[t, x],
  D[v3o1[t, x], 
    t] == (0.001 +10^ttt (v1o1[t, x]^2/(1/k1^2 + v1o1[t, x]^2)) - 
     0.005*(v3o1[t, x]) - 
     0.001*(v3o1[t, x])*((v4o2[t, x] + v4o3[t, x])/2) + 
     0.1*(v5o2[t, x])), 
  D[v4o2[t, x], 
    t] == (-0.001*(v3o1[t, x])*(v4o2[t, x]) + 0.1*(v5o2[t, x]) - 
     0.005*(v4o2[t, x]) + 0.007),
  D[v5o2[t, x], t] == 
   0.001*(v3o1[t, x])*(v4o2[t, x]) - 0.1*(v5o2[t, x]) - 
    0.005*(v5o2[t, x]),
  D[v6o2[t, x], t] == ( 0.1*(v5o2[t, x]) - 0.005*(v6o2[t, x])),
  D[v7o2[t, x], t] == ( 
    0.001 + ((0.1*(v6o2[t, x])^2)/(1 + (v6o2[t, x])^2)) - 
     0.005*(v7o2[t, x])),
  D[v1o2[t, x], t] ==  
   D0 D[v1o2[t, x], x, 
      x] + ((\[Mu]\[Mu]) K1 ((b + v1o2[t, x])/(1 + v1o2[t, x])) v2o2[
       t, x]) - ((\[CapitalGamma]\[CapitalGamma] v1o2[t, x])/(K + 
       v1o2[t, x])) + B - 
    10^zzz(*Q110^xxx/\[Tau]n*)(v7o2[t, x]) v1o2[t, x],
  D[v2o2[t, x], t] == 1/(1 + v1o2[t, x]^2) - v2o2[t, x],
  D[v3o2[t, x], 
    t] == (0.001 +10^ttt (v1o2[t, x]^2/(1/k1^2 + v1o2[t, x]^2)) - 
     0.005*(v3o2[t, x]) - 
     0.001*(v3o2[t, x])*((v4o1[t, x] + v4o3[t, x])/2) + 
     0.1*(v5o1[t, x])), 
  D[v4o1[t, x], 
    t] == (-0.001*(v3o2[t, x])*(v4o1[t, x]) + 0.1*(v5o1[t, x]) - 
     0.005*(v4o1[t, x]) + 0.007),
  D[v5o1[t, x], t] == 
   0.001*(v3o2[t, x])*(v4o1[t, x]) - 0.1*(v5o1[t, x]) - 
    0.005*(v5o1[t, x]),
  D[v6o1[t, x], t] == ( 0.1*(v5o1[t, x]) - 0.005*(v6o1[t, x])),
  D[v7o1[t, x], t] == ( 
    0.001 + ((0.1*(v6o1[t, x])^2)/(1 + (v6o1[t, x])^2)) - 
     0.005*(v7o1[t, x])),
  (**)
  D[v1o3[t, x], t] ==  
   D0 D[v1o3[t, x], x, 
      x] + ((\[Mu]0) K1 ((b + v1o3[t, x])/(1 + v1o3[t, x])) v2o3[t, 
       x]) - ((\[CapitalGamma]\[CapitalGamma] v1o3[t, x])/(K + 
       v1o3[t, x])) + B - 
    10^zzz(*Q110^xxx/\[Tau]n*)(v7o3[t, x]) v1o3[t, x],
  D[v2o3[t, x], t] == 1/(1 + v1o3[t, x]^2) - v2o3[t, x],
  D[v3o3[t, x], 
    t] == (0.001 +10^ttt (v1o3[t, x]^2/(1/k1^2 + v1o3[t, x]^2)) - 
     0.005*(v3o3[t, x]) - 
     0.001*(v3o3[t, x])*((v4o1[t, x] + v4o4[t, x])/2) + 
     0.1*(v5o4[t, x])), 
  D[v4o4[t, x], 
    t] == (-0.001*(v3o3[t, x])*(v4o4[t, x]) + 0.1*(v5o4[t, x]) - 
     0.005*(v4o4[t, x]) + 0.007),
  D[v5o4[t, x], t] == 
   0.001*(v3o3[t, x])*(v4o4[t, x]) - 0.1*(v5o4[t, x]) - 
    0.005*(v5o4[t, x]),
  D[v6o4[t, x], t] == ( 0.1*(v5o4[t, x]) - 0.005*(v6o4[t, x])),
  D[v7o4[t, x], t] == ( 
    0.001 + ((0.1*(v6o4[t, x])^2)/(1 + (v6o4[t, x])^2)) - 
     0.005*(v7o4[t, x])),
  D[v1o4[t, x], t] ==  
   D0 D[v1o4[t, x], x, 
      x] + ((\[Mu]\[Mu]) K1 ((b + v1o4[t, x])/(1 + v1o4[t, x])) v2o4[
       t, x]) - ((\[CapitalGamma]\[CapitalGamma] v1o4[t, x])/(K + 
       v1o4[t, x])) + B - 
    10^zzz(*Q110^xxx/\[Tau]n*)(v7o4[t, x]) v1o4[t, x],
  D[v2o4[t, x], t] == 1/(1 + v1o4[t, x]^2) - v2o4[t, x],
  D[v3o4[t, x], 
    t] == (0.001 +10^ttt (v1o4[t, x]^2/(1/k1^2 + v1o4[t, x]^2)) - 
     0.005*(v3o4[t, x]) - 
     0.001*(v3o4[t, x])*((v4o1[t, x] + v4o3[t, x])/2) + 
     0.1*(v5o3[t, x])), 
  D[v4o3[t, x], 
    t] == (-0.001*(v3o4[t, x])*(v4o3[t, x]) + 0.1*(v5o3[t, x]) - 
     0.005*(v4o3[t, x]) + 0.007),
  D[v5o3[t, x], t] == 
   0.001*(v3o4[t, x])*(v4o3[t, x]) - 0.1*(v5o3[t, x]) - 
    0.005*(v5o3[t, x]),
  D[v6o3[t, x], t] == ( 0.1*(v5o3[t, x]) - 0.005*(v6o3[t, x])),
  D[v7o3[t, x], t] == ( 
    0.001 + ((0.1*(v6o3[t, x])^2)/(1 + (v6o3[t, x])^2)) - 
     0.005*(v7o3[t, x])),
  (**)
  D[v1o5[t, x], t] ==  
   D0 D[v1o5[t, x], x, 
      x] + ((\[Mu]0) K1 ((b + v1o5[t, x])/(1 + v1o5[t, x])) v2o5[t, 
       x]) - ((\[CapitalGamma]\[CapitalGamma] v1o5[t, x])/(K + 
       v1o5[t, x])) + B - 
    10^zzz(*Q110^xxx/\[Tau]n*)(v7o5[t, x]) v1o5[t, x],
  D[v2o5[t, x], t] == 1/(1 + v1o5[t, x]^2) - v2o5[t, x],
  D[v3o5[t, x], 
    t] == (0.001 +10^ttt (v1o5[t, x]^2/(1/k1^2 + v1o5[t, x]^2)) - 
     0.005*(v3o5[t, x]) - 
     0.001*(v3o5[t, x])*((v4o6[t, x] + v4o3[t, x])/2) + 
     0.1*(v5o6[t, x])), 
  D[v4o6[t, x], 
    t] == (-0.001*(v3o5[t, x])*(v4o6[t, x]) + 0.1*(v5o6[t, x]) - 
     0.005*(v4o6[t, x]) + 0.007),
  D[v5o6[t, x], t] == 
   0.001*(v3o5[t, x])*(v4o6[t, x]) - 0.1*(v5o6[t, x]) - 
    0.005*(v5o6[t, x]),
  D[v6 o6[t, x], t] == ( 0.1*(v5o6[t, x]) - 0.005*(v6o6[t, x])),
  D[v7o6[t, x], t] == ( 
    0.001 + ((0.1*(v6o6[t, x])^2)/(1 + (v6o6[t, x])^2)) - 
     0.005*(v7o6[t, x])),
  D[v1o6[t, x], t] ==  
   D0 D[v1o6[t, x], x, 
      x] + ((\[Mu]\[Mu]) K1 ((b + v1o6[t, x])/(1 + v1o6[t, x])) v2o6[
       t, x]) - ((\[CapitalGamma]\[CapitalGamma] v1o6[t, x])/(K + 
       v1o6[t, x])) + B - 
    10^zzz(*Q110^xxx/\[Tau]n*)(v7o6[t, x]) v1o6[t, x],
  D[v2o6[t, x], t] == 1/(1 + v1o6[t, x]^2) - v2o6[t, x],
  D[v3o6[t, x], 
    t] == (0.001 +10^ttt (v1o6[t, x]^2/(1/k1^2 + v1o6[t, x]^2)) - 
     0.005*(v3o6[t, x]) - 
     0.001*(v3o6[t, x])*((v4o5[t, x] + v4o3[t, x])/2) + 
     0.1*(v5o5[t, x])), 
  D[v4o5[t, x], 
    t] == (-0.001*(v3o6[t, x])*(v4o5[t, x]) + 0.1*(v5o5[t, x]) - 
     0.005*(v4o5[t, x]) + 0.007),
  D[v5o5[t, x], t] == 
   0.001*(v3o6[t, x])*(v4o5[t, x]) - 0.1*(v5o5[t, x]) - 
    0.005*(v5o5[t, x]),
  D[v6o5[t, x], t] == ( 0.1*(v5o5[t, x]) - 0.005*(v6o5[t, x])),
  D[v7o5[t, x], t] == ( 
    0.001 + ((0.1*(v6o5[t, x])^2)/(1 + (v6o5[t, x])^2)) - 
     0.005*(v7o5[t, x])),
  v2o1[0, x] == iniope,
  v1o1[0, x] == If[-20 < x < 20, 2 , inical],
  v1o1[t, -200] == inical ,
  v1o1[t, 200] == inical ,
  v3o1[0, x] == 0,
  v4o1[0, x] == 0.1,(**)
  v5o1[0, x] == 0,
  v6o1[0, x] == 0,
  v7o1[0, x] == 0,
  v2o2[0, x] == iniope,
  v1o2[0, x] == If[-20 < x < 20, 2 , inical],
  v1o2[t, -200] == inical ,
  v1o2[t, 200] == inical ,
  v3o2[0, x] == 0,
  v4o2[0, x] == 0.1,(**)
  v5o2[0, x] == 0,
  v6o2[0, x] == 0,
  v7o2[0, x] == 0,
  (**)
  v2o3[0, x] == iniope,
  v1o3[0, x] == If[-20 < x < 20, 2 , inical],
  v1o3[t, -200] == inical ,
  v1o3[t, 200] == inical ,
  v3o3[0, x] == 0,
  v4o3[0, x] == 0.1,(**)
  v5o3[0, x] == 0,
  v6o3[0, x] == 0,
  v7o3[0, x] == 0,
  v2o4[0, x] == iniope,
  v1o4[0, x] == If[-20 < x < 20, 2 , inical],
  v1o4[t, -200] == inical ,
  v1o4[t, 200] == inical ,
  v3o4[0, x] == 0,
  v4o4[0, x] == 0.1,(**)
  v5o4[0, x] == 0,
  v6o4[0, x] == 0,
  v7o4[0, x] == 0,
  (**)
  v2o5[0, x] == iniope,
  v1o5[0, x] == If[-20 < x < 20, 2 , inical],
  v1o5[t, -200] == inical ,
  v1o5[t, 200] == inical ,
  v3o5[0, x] == 0,
  v4o5[0, x] == 0.1,(**)
  v5o5[0, x] == 0,
  v6o5[0, x] == 0,
  v7o5[0, x] == 0,
  v2o6[0, x] == iniope,
  v1o6[0, x] == If[-20 < x < 20, 2 , inical],
  v1o6[t, -200] == inical ,
  v1o6[t, 200] == inical ,
  v3o6[0, x] == 0,
  v4o6[0, x] == 0.1,(**)
  v5o6[0, x] == 0,
  v6o6[0, x] == 0,
  v7o6[0, x] == 0
  },
 {v1o1, v2o1, v3o1, v4o1, v5o1, v6o1, v7o1, v1o2, v2o2, v3o2, v4o2, 
  v5o2, v6o2, v7o2(**), v1o3, v2o3, v3o3, v4o3, v5o3, v6o3, v7o3, 
  v1o4, v2o4, v3o4, v4o4, v5o4, v6o4, v7o4,(**)v1o5, v2o5, v3o5, v4o5,
   v5o5, v6o5, v7o5, v1o6, v2o6, v3o6, v4o6, v5o6, v6o6, v7o6}, {t, 0,
   tmax}, {x, -200, 200},
 Method -> {"MethodOfLines", 
   "SpatialDiscretization" -> {"TensorProductGrid", 
     "MaxPoints" -> 101}}]
Plot[{Evaluate[{v3o1[t, 0]} /. sol], Evaluate[{v3o2[t, 0]} /. sol], 
  Evaluate[{v3o3[t, 0]} /. sol], Evaluate[{v3o4[t, 0]} /. sol], 
  Evaluate[{v3o5[t, 0]} /. sol], Evaluate[{v3o6[t, 0]} /. sol]}, {t, 
  0, tmax}, PlotStyle -> {Red, Blue}, ColorFunction -> "Rainbow", 
 PlotLegends -> Automatic]
```
$\endgroup$
17
  • 2
    $\begingroup$ in theory, with a numerical solver, having 16 equations vs. 32 vs. 64, should not make any difference (except for memory and time). As to why NDSolve can't solve your system when you add more equations, it is hard to know, since you did not show an example or say what error or what happens exactly. does it hang? what does it do? Without an example, hard to say. $\endgroup$
    – Nasser
    Commented Feb 28, 2021 at 9:28
  • $\begingroup$ @Nasser you are right, I added the code. $\endgroup$
    – confused
    Commented Feb 28, 2021 at 10:15
  • 2
    $\begingroup$ @AlexTrounev That's a PDE, isn't it? An ODE has the form u'[t] == F[t, u[t]....I think we should be strict in our definitions when dealing with a computer program, or at least as strict as the program. $\endgroup$
    – Michael E2
    Commented Feb 28, 2021 at 15:00
  • 1
    $\begingroup$ @confused Quitting the kernel is known to fix inexplicable problems. If it works, it should be because of a forgotten definition that wasn't cleared. $\endgroup$
    – Michael E2
    Commented Feb 28, 2021 at 15:01
  • 2
    $\begingroup$ Quitting the kernel could also reset a bug. Like rebooting the computer. :) $\endgroup$
    – Michael E2
    Commented Feb 28, 2021 at 15:15

1 Answer 1

4
$\begingroup$

After some debugging this code works on v12.0 and 12.2 as well

Clear["Global`*"]
(*SET UP*)
tmax = 2500;
(*INITIAL CONDITIONS*)

inical = SetPrecision[0.200165185171545551`5., 5];
iniope = 1/(1 + inical^2);
\[Mu]\[Mu] = SetPrecision[0.5500000, 5];
\[Mu]0 = SetPrecision[0.5155877999949406`, 5];
(*PARAMETERS*)
Dc = 20;
kflux = 8.1;
l = 20;(*--------*)b = 0.111 // Rationalize;
k1 = 0.7;
\[Gamma] = 2.0;
k\[Gamma] = 0.1;
\[Beta] =(*0*)0.02;
\[Tau]n = 2.0;
k2 = 0.7;
K1 = ((kflux \[Tau]n)/k1);
\[CapitalGamma]\[CapitalGamma] = ((\[Tau]n \[Gamma] 0.9)/k1);
\[CapitalGamma] = ((\[Tau]n \[Gamma])/k1);
T = (\[Gamma]*\[Tau]n)/k1;
K = (k\[Gamma]/k1);
K2 = k2/k1;
B = \[Beta]/k1;
D0 = Dc \[Tau]n/l^2;
Dp0 = Dp \[Tau]n/l^2;
Dp = 300;
\[Mu]1 = 0.433;
kp = 0.02;
(*q=-3.4;*)
xxx = -0.1;
qqq = -0.4;
Q1 = 0.09;
Q2 = 0.5;
ttt = 0.6;
zzz = -1.2; var = {v1o1, v2o1, v3o1, v4o1, v5o1, v6o1, v7o1, v1o2, 
  v2o2, v3o2, v4o2, v5o2, v6o2, v7o2, v1o3, v2o3, v3o3, v4o3, v5o3, 
  v6o3, v7o3, v1o4, v2o4, v3o4, v4o4, v5o4, v6o4, v7o4, v1o5, v2o5, 
  v3o5, v4o5, v5o5, v6o5, v7o5, v1o6, v2o6, v3o6, v4o6, v5o6, v6o6, 
  v7o6};
eq = {D[v1o1[t, x], t] == 
    D0 D[v1o1[t, x], x, 
       x] + ((\[Mu]0) K1 ((b + v1o1[t, x])/(1 + v1o1[t, x])) v2o1[t, 
        x]) - ((\[CapitalGamma]\[CapitalGamma] v1o1[t, x])/(K + 
         v1o1[t, x])) + B - 10^zzz (v7o1[t, x]) v1o1[t, x], 
   D[v2o1[t, x], t] == 1/(1 + v1o1[t, x]^2) - v2o1[t, x], 
   D[v3o1[t, x], 
     t] == (0.001 + 10^ttt (v1o1[t, x]^2/(1/k1^2 + v1o1[t, x]^2)) - 
      0.005*(v3o1[t, x]) - 
      0.001*(v3o1[t, x])*((v4o2[t, x] + v4o3[t, x])/2) + 
      0.1*(v5o2[t, x])), 
   D[v4o2[t, x], 
     t] == (-0.001*(v3o1[t, x])*(v4o2[t, x]) + 0.1*(v5o2[t, x]) - 
      0.005*(v4o2[t, x]) + 0.007), 
   D[v5o2[t, x], t] == 
    0.001*(v3o1[t, x])*(v4o2[t, x]) - 0.1*(v5o2[t, x]) - 
     0.005*(v5o2[t, x]), 
   D[v6o2[t, x], t] == (0.1*(v5o2[t, x]) - 0.005*(v6o2[t, x])), 
   D[v7o2[t, x], 
     t] == (0.001 + ((0.1*(v6o2[t, x])^2)/(1 + (v6o2[t, x])^2)) - 
      0.005*(v7o2[t, x])), 
   D[v1o2[t, x], t] == 
    D0 D[v1o2[t, x], x, 
       x] + ((\[Mu]\[Mu]) K1 ((b + v1o2[t, x])/(1 + v1o2[t, x])) v2o2[
        t, x]) - ((\[CapitalGamma]\[CapitalGamma] v1o2[t, x])/(K + 
         v1o2[t, x])) + B - 10^zzz (v7o2[t, x]) v1o2[t, x], 
   D[v2o2[t, x], t] == 1/(1 + v1o2[t, x]^2) - v2o2[t, x], 
   D[v3o2[t, x], 
     t] == (0.001 + 10^ttt (v1o2[t, x]^2/(1/k1^2 + v1o2[t, x]^2)) - 
      0.005*(v3o2[t, x]) - 
      0.001*(v3o2[t, x])*((v4o1[t, x] + v4o3[t, x])/2) + 
      0.1*(v5o1[t, x])), 
   D[v4o1[t, x], 
     t] == (-0.001*(v3o2[t, x])*(v4o1[t, x]) + 0.1*(v5o1[t, x]) - 
      0.005*(v4o1[t, x]) + 0.007), 
   D[v5o1[t, x], t] == 
    0.001*(v3o2[t, x])*(v4o1[t, x]) - 0.1*(v5o1[t, x]) - 
     0.005*(v5o1[t, x]), 
   D[v6o1[t, x], t] == (0.1*(v5o1[t, x]) - 0.005*(v6o1[t, x])), 
   D[v7o1[t, x], 
     t] == (0.001 + ((0.1*(v6o1[t, x])^2)/(1 + (v6o1[t, x])^2)) - 
      0.005*(v7o1[t, x])),(**)
   D[v1o3[t, x], t] == 
    D0 D[v1o3[t, x], x, 
       x] + ((\[Mu]0) K1 ((b + v1o3[t, x])/(1 + v1o3[t, x])) v2o3[t, 
        x]) - ((\[CapitalGamma]\[CapitalGamma] v1o3[t, x])/(K + 
         v1o3[t, x])) + B - 10^zzz (v7o3[t, x]) v1o3[t, x], 
   D[v2o3[t, x], t] == 1/(1 + v1o3[t, x]^2) - v2o3[t, x], 
   D[v3o3[t, x], 
     t] == (0.001 + 10^ttt (v1o3[t, x]^2/(1/k1^2 + v1o3[t, x]^2)) - 
      0.005*(v3o3[t, x]) - 
      0.001*(v3o3[t, x])*((v4o1[t, x] + v4o4[t, x])/2) + 
      0.1*(v5o4[t, x])), 
   D[v4o4[t, x], 
     t] == (-0.001*(v3o3[t, x])*(v4o4[t, x]) + 0.1*(v5o4[t, x]) - 
      0.005*(v4o4[t, x]) + 0.007), 
   D[v5o4[t, x], t] == 
    0.001*(v3o3[t, x])*(v4o4[t, x]) - 0.1*(v5o4[t, x]) - 
     0.005*(v5o4[t, x]), 
   D[v6o4[t, x], t] == (0.1*(v5o4[t, x]) - 0.005*(v6o4[t, x])), 
   D[v7o4[t, x], 
     t] == (0.001 + ((0.1*(v6o4[t, x])^2)/(1 + (v6o4[t, x])^2)) - 
      0.005*(v7o4[t, x])), 
   D[v1o4[t, x], t] == 
    D0 D[v1o4[t, x], x, 
       x] + ((\[Mu]\[Mu]) K1 ((b + v1o4[t, x])/(1 + v1o4[t, x])) v2o4[
        t, x]) - ((\[CapitalGamma]\[CapitalGamma] v1o4[t, x])/(K + 
         v1o4[t, x])) + B - 10^zzz (v7o4[t, x]) v1o4[t, x], 
   D[v2o4[t, x], t] == 1/(1 + v1o4[t, x]^2) - v2o4[t, x], 
   D[v3o4[t, x], 
     t] == (0.001 + 10^ttt (v1o4[t, x]^2/(1/k1^2 + v1o4[t, x]^2)) - 
      0.005*(v3o4[t, x]) - 
      0.001*(v3o4[t, x])*((v4o1[t, x] + v4o3[t, x])/2) + 
      0.1*(v5o3[t, x])), 
   D[v4o3[t, x], 
     t] == (-0.001*(v3o4[t, x])*(v4o3[t, x]) + 0.1*(v5o3[t, x]) - 
      0.005*(v4o3[t, x]) + 0.007), 
   D[v5o3[t, x], t] == 
    0.001*(v3o4[t, x])*(v4o3[t, x]) - 0.1*(v5o3[t, x]) - 
     0.005*(v5o3[t, x]), 
   D[v6o3[t, x], t] == (0.1*(v5o3[t, x]) - 0.005*(v6o3[t, x])), 
   D[v7o3[t, x], 
     t] == (0.001 + ((0.1*(v6o3[t, x])^2)/(1 + (v6o3[t, x])^2)) - 
      0.005*(v7o3[t, x])),(**)
   D[v1o5[t, x], t] == 
    D0 D[v1o5[t, x], x, 
       x] + ((\[Mu]0) K1 ((b + v1o5[t, x])/(1 + v1o5[t, x])) v2o5[t, 
        x]) - ((\[CapitalGamma]\[CapitalGamma] v1o5[t, x])/(K + 
         v1o5[t, x])) + B - 10^zzz (v7o5[t, x]) v1o5[t, x], 
   D[v2o5[t, x], t] == 1/(1 + v1o5[t, x]^2) - v2o5[t, x], 
   D[v3o5[t, x], 
     t] == (0.001 + 10^ttt (v1o5[t, x]^2/(1/k1^2 + v1o5[t, x]^2)) - 
      0.005*(v3o5[t, x]) - 
      0.001*(v3o5[t, x])*((v4o6[t, x] + v4o3[t, x])/2) + 
      0.1*(v5o6[t, x])), 
   D[v4o6[t, x], 
     t] == (-0.001*(v3o5[t, x])*(v4o6[t, x]) + 0.1*(v5o6[t, x]) - 
      0.005*(v4o6[t, x]) + 0.007), 
   D[v5o6[t, x], t] == 
    0.001*(v3o5[t, x])*(v4o6[t, x]) - 0.1*(v5o6[t, x]) - 
     0.005*(v5o6[t, x]), 
   D[v6o6[t, x], t] == (0.1*(v5o6[t, x]) - 0.005*(v6o6[t, x])), 
   D[v7o6[t, x], 
     t] == (0.001 + ((0.1*(v6o6[t, x])^2)/(1 + (v6o6[t, x])^2)) - 
      0.005*(v7o6[t, x])), 
   D[v1o6[t, x], t] == 
    D0 D[v1o6[t, x], x, 
       x] + ((\[Mu]\[Mu]) K1 ((b + v1o6[t, x])/(1 + v1o6[t, x])) v2o6[
        t, x]) - ((\[CapitalGamma]\[CapitalGamma] v1o6[t, x])/(K + 
         v1o6[t, x])) + B - 10^zzz (v7o6[t, x]) v1o6[t, x], 
   D[v2o6[t, x], t] == 1/(1 + v1o6[t, x]^2) - v2o6[t, x], 
   D[v3o6[t, x], 
     t] == (0.001 + 10^ttt (v1o6[t, x]^2/(1/k1^2 + v1o6[t, x]^2)) - 
      0.005*(v3o6[t, x]) - 
      0.001*(v3o6[t, x])*((v4o5[t, x] + v4o3[t, x])/2) + 
      0.1*(v5o5[t, x])), 
   D[v4o5[t, x], 
     t] == (-0.001*(v3o6[t, x])*(v4o5[t, x]) + 0.1*(v5o5[t, x]) - 
      0.005*(v4o5[t, x]) + 0.007), 
   D[v5o5[t, x], t] == 
    0.001*(v3o6[t, x])*(v4o5[t, x]) - 0.1*(v5o5[t, x]) - 
     0.005*(v5o5[t, x]), 
   D[v6o5[t, x], t] == (0.1*(v5o5[t, x]) - 0.005*(v6o5[t, x])), 
   D[v7o5[t, x], 
     t] == (0.001 + ((0.1*(v6o5[t, x])^2)/(1 + (v6o5[t, x])^2)) - 
      0.005*(v7o5[t, x]))};
icbc = {v2o1[0, x] == iniope, 
   v1o1[0, x] == If[-20 < x < 20, 2, inical], v1o1[t, -200] == inical,
    v1o1[t, 200] == inical, v3o1[0, x] == 0, v4o1[0, x] == 0.1,(**)
   v5o1[0, x] == 0, v6o1[0, x] == 0, v7o1[0, x] == 0, 
   v2o2[0, x] == iniope, v1o2[0, x] == If[-20 < x < 20, 2, inical], 
   v1o2[t, -200] == inical, v1o2[t, 200] == inical, v3o2[0, x] == 0, 
   v4o2[0, x] == 0.1,(**)v5o2[0, x] == 0, v6o2[0, x] == 0, 
   v7o2[0, x] == 0,(**)v2o3[0, x] == iniope, 
   v1o3[0, x] == If[-20 < x < 20, 2, inical], v1o3[t, -200] == inical,
    v1o3[t, 200] == inical, v3o3[0, x] == 0, v4o3[0, x] == 0.1,(**)
   v5o3[0, x] == 0, v6o3[0, x] == 0, v7o3[0, x] == 0, 
   v2o4[0, x] == iniope, v1o4[0, x] == If[-20 < x < 20, 2, inical], 
   v1o4[t, -200] == inical, v1o4[t, 200] == inical, v3o4[0, x] == 0, 
   v4o4[0, x] == 0.1,(**)v5o4[0, x] == 0, v6o4[0, x] == 0, 
   v7o4[0, x] == 0,(**)v2o5[0, x] == iniope, 
   v1o5[0, x] == If[-20 < x < 20, 2, inical], v1o5[t, -200] == inical,
    v1o5[t, 200] == inical, v3o5[0, x] == 0, v4o5[0, x] == 0.1,(**)
   v5o5[0, x] == 0, v6o5[0, x] == 0, v7o5[0, x] == 0, 
   v2o6[0, x] == iniope, v1o6[0, x] == If[-20 < x < 20, 2, inical], 
   v1o6[t, -200] == inical, v1o6[t, 200] == inical, v3o6[0, x] == 0, 
   v4o6[0, x] == 0.1,(**)v5o6[0, x] == 0, v6o6[0, x] == 0, 
   v7o6[0, x] == 0};
sol = NDSolve[{eq, icbc}, var, {t, 0, tmax}, {x, -200, 200}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MaxPoints" -> 101}}];

Visualization of numerical solution

Table[DensityPlot[
  Evaluate[var[[i]][t, x] /. sol[[1]]], {t, 0, tmax}, {x, -200, 200}, 
  PlotLabel -> var[[i]], ColorFunction -> "Rainbow", 
  PlotLegends -> Automatic, PlotRange -> All], {i, Length[var]}]

Figure 1

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.