0
$\begingroup$

I'm trying to numerically solve the following differential equations eqs for plotting the bifurcation diagram. However, under special values of the parameter u, I get a problem when calculating the differential equations with NDSolve. The corresponding codes are shown as follows

Initial definitions

\[Alpha] = 0.005;\[Beta] = 0.2;\[Gamma] = 10;kt = 100;\[Xi]b = 0.82;
\[Beta]L = {1.8751040687119611791`8., 4.6940911330190088767`8.,7.8547574382376149826`8., 10.9955407348773432322`8.,14.1371683910464717872`8., 17.2787595320882612395`8.,20.4203522521036040901`8.};
Subscript[\[Phi], m_][x] :=Cos[\[Beta]L[[m]]*x] - Cosh[\[Beta]L[[m]]*x] + ((Sin[\[Beta]L[[m]]] - Sinh[\[Beta]L[[m]]])/(Cos[\[Beta]L[[m]]] + Cosh[\[Beta]L[[m]]]))*(Sin[\[Beta]L[[m]]*x] - Sinh[\[Beta]L[[m]]*x]);
nig[fx_] :=N[Round[Chop[NIntegrate[fx, {x, 0, 1}, AccuracyGoal -> 10, WorkingPrecision -> 40, MaxRecursion -> 20], 10^-8], 10^-6], 6];

Differential equations

mm = 4;
eqs = Table[\!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\((\[Alpha]*nig[\(\*SubscriptBox[\(\[Phi]\), \(i\)]\)[x]*D[\(\*SubscriptBox[\(\[Phi]\), \(j\)]\)[x], {x, 4}]]*\(\*SubscriptBox[\(q\), \(j\)]'\)[t])\)\) 
      +\!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\((nig[\(\*SubscriptBox[\(\[Phi]\), \(i\)]\)[x]*D[\(\*SubscriptBox[\(\[Phi]\), \(j\)]\)[x], {x, 4}]]*\(\*SubscriptBox[\(q\), \(j\)]\)[t])\)\)
      +\!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\((nig[\(\*SubscriptBox[\(\[Phi]\), \(i\)]\)[x]*D[\(\*SubscriptBox[\(\[Phi]\), \(j\)]\)[x], {x, 2}]]*\*SuperscriptBox[\(u\), \(2\)]*\(\*SubscriptBox[\(q\), \(j\)]\)[t])\)\)
      -\!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\((nig[\(\*SubscriptBox[\(\[Phi]\), \(i\)]\)[x]*D[\(\*SubscriptBox[\(\[Phi]\), \(j\)]\)[x], {x, 2}]*\[Gamma]*\((1 - x)\)]*\(\*SubscriptBox[\(q\), \(j\)]\)[t])\)\)
      +\!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\((nig[\(\*SubscriptBox[\(\[Phi]\), \(i\)]\)[x]*D[\(\*SubscriptBox[\(\[Phi]\), \(j\)]\)[x], x]]*2*\*SqrtBox[\(\[Beta]\)]*u*\(\*SubscriptBox[\(q\), \(j\)]'\)[t])\)\)
      +\!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\((nig[\(\*SubscriptBox[\(\[Phi]\), \(i\)]\)[x]*D[\(\*SubscriptBox[\(\[Phi]\), \(j\)]\)[x], x]]*\[Gamma]*\(\*SubscriptBox[\(q\), \(j\)]\)[t])\)\)
      +\!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\(\*UnderoverscriptBox[\(\[Sum]\), \(k = 1\), \(mm\)]\(\*UnderoverscriptBox[\(\[Sum]\), \(l = 1\), \(mm\)]\((\((\(\*SubscriptBox[\(\[Phi]\), \(i\)]\)[x] /. x -> \[Xi]b)\)*\((\(\*SubscriptBox[\(\[Phi]\), \(j\)]\)[x] /. x -> \[Xi]b)\)*\((\(\*SubscriptBox[\(\[Phi]\), \(k\)]\)[x] /. x -> \[Xi]b)\)*\((\(\*SubscriptBox[\(\[Phi]\), \(l\)]\)[x] /. x -> \[Xi]b)\)*kt*\(\*SubscriptBox[\(q\), \(j\)]\)[t]*\(\*SubscriptBox[\(q\), \(k\)]\)[t]*\(\*SubscriptBox[\(q\), \(l\)]\)[t])\)\)\)\)
      +\!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\((nig[\(\*SubscriptBox[\(\[Phi]\), \(i\)]\)[x]*\(\*SubscriptBox[\(\[Phi]\), \(j\)]\)[x]]*\(\*SubscriptBox[\(q\), \(j\)]''\)[t])\)\) == 0, {i, mm}];
EQs = Chop[Expand[Join[eqs]], 10^-8]

Initial conditions

{q101, q102, q201, q202, q301, q302, q401, q402} = {0.001, 0, 0, 0, 0,
0, 0, 0};
ICs = {Subscript[q, 1][0] == q101, Subscript[q, 1]'[0] == q102, Subscript[q, 2][0] == q201, Subscript[q, 2]'[0] == q202, Subscript[q, 3][0] == q301, Subscript[q, 3]'[0] == q302, Subscript[q, 4][0] == q401, Subscript[q, 4]'[0] == q402};


u = 8.8;
s0 = NDSolve[Join[EQs, ICs], {Subscript[q, 1], Subscript[q, 2], Subscript[q, 3], Subscript[q, 4]}, {t, 0, 1000}, Method -> "ImplicitRungeKutta"]

The problem is that when the parameter value u is lower than 8.8, like u=8.7 or even u=6.5, the above-mentioned codes are efficient for the calculation in the interval $u\in\left(6,8.7\right)$. But when the parameter u beyonds 8.7, such as u=8.8, 8.9, and 9, etc, I get the following warning:

NDSolve::mxst: Maximum number of 216267 steps reached at the point t == 83.2334874133867`

Concerning this warning, I use

  • MaxSteps -> 1000000
  • WorkingPrecision -> 40

But Mathematica still gives me the warning. When selecting MaxSteps -> Infinity, Mathematica will keeps running for a long time. And then, I stops it.


Question

I would like to numerically solve the differential equations eqs with $u\in\left(6,9.5\right)$. So I am wondering:

  • Is there a specific method for numerically solving the differential equations? It may be not the "ImplicitRungeKutta" or even the "NDSolve".
  • I succeeded in calculating these equations when choosing mm=2, where mm is the number of 2-order differential equations considered. I suspect that the appearance of the warning NDSolve::mxst: Maximum number of $\times\times$ steps reached at the point t == $\times\times$ may have something to do with the value of mm, but I am not sure. Is there any other way I can do that when mm>4?

Thanks for help and insights in advance!


Supplement

I have noticed MarcoB's suggestions. Some related replies could be seen in the comment section. Here, I add some changes to the subscript format of the code.

Initial definitions

\[Phi][m_, x_] := Cos[\[Beta]L[[m]]*x] - Cosh[\[Beta]L[[m]]*x] + ((Sin[\[Beta]L[[m]]] - Sinh[\[Beta]L[[m]]])/(Cos[\[Beta]L[[m]]] + Cosh[\[Beta]L[[m]]]))*(Sin[\[Beta]L[[m]]*x] - Sinh[\[Beta]L[[m]]*x]);

Differential equations

eqs = Table[\!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\((\
[Alpha]*nig[\[Phi][i, x]*D[\[Phi][j, x], {x, 4}]]*D[q[j, t], t])\)\) 
+ \!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\((nig[\[Phi][i, 
x]*D[\[Phi][j, x], {x, 4}]]*q[j, t])\)\) 
+ \!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\((nig[\[Phi][i, 
x]*D[\[Phi][j, x], {x, 2}]]*\*SuperscriptBox[\(u\), \(2\)]*q[j, t])\)\) 
- \!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\((nig[\[Phi][i, 
x]*D[\[Phi][j, x], {x, 2}]*\[Gamma]*\((1 - x)\)]*q[j, t])\)\) 
+ \!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\((nig[\[Phi][i, 
x]*D[\[Phi][j, x], x]]*2*\*SqrtBox[\(\[Beta]\)]*u*D[q[j, t], t])\)\) 
+ \!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\((nig[\[Phi][i, 
x]*D[\[Phi][j, x], x]]*\[Gamma]*q[j, t])\)\) 
+ \!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\
(\*UnderoverscriptBox[\(\[Sum]\), \(k = 1\), \(mm\)]\(\*UnderoverscriptBox[\
(\[Sum]\), \(l = 1\), \(mm\)]\((\((\[Phi][i, x] /. x -> \[Xi]b)\)*\((\[Phi]
[j, x] /. x -> \[Xi]b)\)*\((\[Phi][k, x] /. x -> \[Xi]b)\)*\((\[Phi][l, x] /. 
x -> \[Xi]b)\)*kt*q[j, t]*q[k, t]*q[l, t])\)\)\)\) 
+ \!\(\*UnderoverscriptBox[\(\[Sum]\), \(j = 1\), \(mm\)]\((nig[\[Phi][i, 
x]*\[Phi][j, x]]*D[q[j, t], {t, 2}])\)\) == 0, {i, mm}];

EQs = EQs /. {q[1, t] -> q1[t], q[2, t] -> q2[t], q[3, t] -> q3[t], q[4, t] -> q4[t], \!\(\*SuperscriptBox[\(q\), TagBox[RowBox[{"(", RowBox[{"0", ",", "1"}], ")"}],Derivative],MultilineFunction->None]\)[1, t] -> q1'[t], \!(\*SuperscriptBox[\(q\), TagBox[RowBox[{"(", RowBox[{"0", ",", "1"}], ")"}],Derivative],MultilineFunction->None]\)[2, t] -> q2'[t], \!\(\*SuperscriptBox[\(q\), TagBox[RowBox[{"(", RowBox[{"0", ",", "1"}], ")"}],Derivative],MultilineFunction->None]\)[3, t] -> q3'[t], \!\(\*SuperscriptBox[\(q\), TagBox[RowBox[{"(", RowBox[{"0", ",", "1"}], ")"}],Derivative],MultilineFunction->None]\)[4, t] -> q4'[t], \!\(\*SuperscriptBox[\(q\), TagBox[RowBox[{"(", RowBox[{"0", ",", "2"}], ")"}],Derivative],MultilineFunction->None]\)[1, t] -> q1''[t], \!\(\*SuperscriptBox[\(q\), TagBox[RowBox[{"(", RowBox[{"0", ",", "2"}], ")"}],Derivative],MultilineFunction->None]\)[2, t] -> q2''[t], \!\(\*SuperscriptBox[\(q\), TagBox[RowBox[{"(", RowBox[{"0", ",", "2"}], ")"}],Derivative],MultilineFunction->None]\)[3, t] -> q3''[t], \!\(\*SuperscriptBox[\(q\), TagBox[RowBox[{"(", RowBox[{"0", ",", "2"}], ")"}],Derivative],MultilineFunction->None]\)[4, t] -> q4''[t]}

Initial conditions

q1[0] == 0.001, q1'[0] == 0, q2[0] == 0, q2'[0] == 0, q3[0] == 0, q3'[0] == 0, q4[0] == 0, q4'[0] == 0
$\endgroup$
  • 1
    $\begingroup$ Did you notice that your initial conditions are defined as ICs, but then in your NDSolve call you use an undefined ICS1? Is that a typo? Also, your betaL values are defined with 8 digits of precision; it does not make sense to ask for a WorkingPrecision of 40 later on. Have you tried removing all the WorkingPrecision, and explicit Method calls? What happens if you let NDSolve choose the method automatically? Also, Subscripts are mostly a formatting convenience, but they are best avoided in calculation, where they might give rise to sneaky problems. $\endgroup$ – MarcoB Dec 23 '19 at 16:08
  • $\begingroup$ @MarcoB thanks. Firstly, I would change the ICS1 to ICS later. It was right in my Mathematica codes before. Secondly, when the "WorkingPrecision->40" is neglected, the warning "the numerical integral converges too slowly" will happends. Changing the Subscripts and removing the WorkingPrecision or the explicit method, the warning "NDSolve::mxst: Maximum number of ×××× steps reached at the point t == ××" still has. Could it be something else? $\endgroup$ – kuzb Dec 24 '19 at 2:50
  • $\begingroup$ @MarcoB I have tried to numerically solve the differential equations with mm=2 using the ode45 of Matlab. At u=8.8, 8.9, and 9, Matlab could easily solve it. The possibility of the problem may comes from the Method. ode45-NDSolve $\endgroup$ – kuzb Dec 24 '19 at 13:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.