1
$\begingroup$

There might be some relevant questions and answers, but I couldn't complete my code by myself just by reading them. Sorry that I am very new to Mathematica and my question could be a bit rough and not well-organized.

I am trying to find the parameter sets of my differential equations that satisfy some conditions I want to impose.

For a much simpler example, say these are my equations of interest,

s = NDSolve[{a y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 10}];
ss = NDSolve[{b y''[x] + c Sin[y[x]] y[x] == 0, y[0] == 1, y'[0] == 0}, y, {x, 0,10}];

and I want to find the sets (a,b,c) that satisfy

y[6]/.s==y[6]/.ss

For instance,

a = -1;
b = 1.4882301;
c = 3;

this set roughly satisfies the condition but is presumably not the only solution. So what I want to do is, if I simply fix c=3, by tracking a and b, I will be able to plot a-b that the condition is satisfied. Likewise, I want to find all the sets (maybe sequences of sets) under certain conditions for the parameters (fixed a, fixed b,...).

It may take a huge amount of time, but basically, I can do this with my hand by finding one solution and varying the parameters infinitesimally. The point is, I want this to be done as automatically as possible. It might require some steps to be there, but my final aim is, put c=3 and press enter are all that I need to do once the code is made.

I feel like in some way I should use Module (and maybe FindRoot?), but not sure about how to. Can anyone give me some help on how to achieve this? Again, I am very poor at this software yet and even just giving me the right direction will be appreciated.

Thank you


Thanks to @BobHanlon I could catch somehow how to approach this problem. However, if I apply this method to a more complicated version, it doesn't work. Here are the actual functions that I want to work on

eqold[\[CapitalSigma]_] = {Z'[s] == -Sin[\[Psi][s]], \[Psi]''[
    s] == -\[Psi]'[s]/X[s] Cos[\[Psi][s]] + 
    Cos[\[Psi][s]] Sin[\[Psi][s]]/X[s]^2 + \[Gamma][s]/
      X[s] Sin[\[Psi][s]] + P X[s]/2 Cos[\[Psi][s]], \[Gamma]'[
    s] == (\[Psi]'[s] - C0)^2/2 - Sin[\[Psi][s]]^2/(2 X[s]^2) + 
    P X[s] Sin[\[Psi][s]] + \[CapitalSigma], X'[s] == Cos[\[Psi][s]]};
help[eq_] := (Subtract @@ eq // Together // Numerator) == 0
C0 = 0; s1 = 2 \[Pi];
nobsoln = 
 ParametricNDSolve[{eqold[\[CapitalSigma]0][[{1, 4}]], 
   help /@ eqold[\[CapitalSigma]0][[{2, 3}]], \[Psi][0] == 
    0, \[Psi]'[0] == a, X[0] == 0, \[Gamma][0] == 0, 
   Z[0] == 0}, {X, \[Psi], Z, \[Gamma]}, {s, 0, s1}, {a, 
   P, \[CapitalSigma]0}, SolveDelayed -> True]; sobsoln = 
 ParametricNDSolve[{eqold[\[CapitalSigma]0][[{1, 4}]], 
   help /@ eqold[\[CapitalSigma]0][[{2, 3}]], \[Psi][
     s1] == \[Pi], \[Psi]'[s1] == a, X[s1] == 0, \[Gamma][s1] == 0, 
   Z[s1] == 0}, {X, \[Psi], Z, \[Gamma]}, {s, 0, s1}, {a, 
   P, \[CapitalSigma]0}, SolveDelayed -> True];

(Sorry for this strange display of the greek letters, I simply copied and pasted the code and am not sure how to fix this).

The matching conditions I want to impose are

X[3] /. nobsoln == X[3] /. sobsoln
\[Psi][3] /. nobsoln == \[Psi][3] /. sobsoln

Or actually, I want nobsoln and sobsoln to be identical except for s=0 and s=s1, where they are singular.

I know that

P = 1.37889;
\[CapitalSigma]0 = -1.1 P^(2/3);
a = 0.0895105;

is one of the solutions.

Then I defined the "data" like you did (this was typically aimed to check the above solution)

data = Table[{\[CapitalSigma]0 /. 
      FindRoot[(X[a, 1.37889, \[CapitalSigma]0][3] /. 
          nobsoln) == (X[a, 1.37889, \[CapitalSigma]0][3] /. 
          sobsoln), {\[CapitalSigma]0, -1.36}], a}, {a, 
     Range[0.07, 0.09, 10^(-4)]}] // Quiet;

The issue pops up here. If I run this, the "running" continues for some seconds and suddenly just stops and resets things like I "quit kernel". This kind of error happened occasionally for some other codes, but they were resolved just by rerunning, but this one never works however I try.

Can someone figure out why this happens either in general or typically for this code? and give me some useful advice on how to fix it?

$\endgroup$
7
  • $\begingroup$ The problem is singular in a case of actual systems with initial points X[0] == 0 and X[s1] == 0 for nobsoln and sobsoln consequently. This singularity crashes kernel without message. $\endgroup$ Commented Nov 22, 2021 at 4:05
  • $\begingroup$ @AlexTrounev Thanks for pointing out the origin of this issue. Can you also give me some advice on how to fix this? $\endgroup$
    – physgj
    Commented Nov 22, 2021 at 11:34
  • $\begingroup$ I don't understand your problem as you want here "I want nobsoln and sobsoln to be identical except for s=0 and s=s1, where they are singular". Then what do you mean here X[3] /. nobsoln == X[3] /. sobsoln; \[Psi][3] /. nobsoln == \[Psi][3] /. sobsoln? $\endgroup$ Commented Nov 22, 2021 at 15:48
  • $\begingroup$ @AlexTrounev Actually, s=3 is what I picked arbitrarily. Basically, for the desired parameter sets, nobsoln and sobsoln are identical except for s=0 and s=s1. So I am checking if they are equal by examining an arbitrary point of s. $\endgroup$
    – physgj
    Commented Nov 22, 2021 at 16:23
  • $\begingroup$ This code may help a bit in understanding what I mean by p = 1.37889; \[Sigma]0 = -1.1 p^(2/3); C0 = 0; aa = -0.08846200000000001 p/\[Sigma]0; s1 = 2 \[Pi]; Show[Plot[ Evaluate[X[aa, p, \[Sigma]0][s] /. nobsoln], {s, 0, 2 \[Pi]}, PlotStyle -> Blue], Plot[Evaluate[\[Psi][aa, p, \[Sigma]0][s] /. nobsoln], {s, 0, 2 \[Pi]}, PlotStyle -> Blue], Plot[Evaluate[X[aa, p, \[Sigma]0][s] /. sobsoln], {s, 0, 2 \[Pi]}, PlotStyle -> {Red, Dashed}], Plot[Evaluate[\[Psi][aa, p, \[Sigma]0][s] /. sobsoln], {s, 0, 2 \[Pi]}, PlotStyle -> {Red, Dashed}], PlotRange -> Full] $\endgroup$
    – physgj
    Commented Nov 22, 2021 at 16:24

1 Answer 1

1
$\begingroup$
Clear["Global`*"]

$Version

(* "12.3.1 for Mac OS X x86 (64-bit) (June 19, 2021)" *)

eqns1 = {a y'[x] == y[x] Cos[x + y[x]], y[0] == 1};

Use ParametricNDSolve

s = ParametricNDSolve[eqns1, y, {x, 0, 10}, {a}];

plt1 = Plot[
   Evaluate@
    Table[
     Tooltip[y[a][x] /. s, a],
     {a, {-2, -1, 1, 2}}],
   {x, 0, 10},
   AxesLabel -> (Style[#, 12, Bold] & /@ {x, y}),
   PlotLegends -> LineLegend[{-2, -1, 1, 2},
     LegendLabel -> Style[a, 12, Bold]]];

Similarly,

eqns2[b_, c_] = {b y''[x] + c Sin[y[x]] y[x] == 0, y[0] == 1, y'[0] == 0};

ss = ParametricNDSolve[eqns2[b, c], y, {x, 0, 10}, {b, c}];

With[{c = 3},
  plt2 = Plot[
    Evaluate@
     Table[
      Tooltip[y[b, c][x] /. ss, b],
      {b, {-2, -1, 1, 2}}],
    {x, 0, 10},
    PlotStyle -> Dashed,
    AxesLabel -> (Style[#, 12, Bold] & /@ {x, y}),
    PlotLegends -> LineLegend[{-2, -1, 1, 2},
      LegendLabel -> Style[b, 12, Bold]]]];

Combining the plots,

Show[plt1, plt2, PlotRange -> All, AspectRatio -> 1]

enter image description here

To find the intersections for a fixed value of c (e.g., c == 3) use FindRoot. Note that good initial estimates for a parameter must be given to FindRoot

data =
  Table[
    Table[
     {a, b /. FindRoot[(y[a][6] /. s) == (y[b, 3][6] /. ss), {b, est}]},
     {a, Range[-2, 2, 0.1] /. {0. -> Nothing}}],
    {est, 0.1, 2, 0.1}] // Quiet;

Plotting the data,

ListLinePlot[data, PlotStyle -> {Automatic, Dashed},
 AxesLabel -> (Style[#, 12, Bold] & /@ {a, b})]

enter image description here

For a == -1 the found values (i.e., may not be all-inclusive) of b are

Last /@ (
  Mean /@
   GatherBy[
    Cases[data, {-1., _}, Infinity],
    Round[#[[2]], 0.000001] &])

(* {0.108524, 0.232989, 0.314367, 0.814686, 1.48823} *)
$\endgroup$
3
  • $\begingroup$ Thank you very much for your help. I think this method makes sense. I really appreciate it. However, if I try to apply the same strategy to my actual functions, it keeps failing to run the "data" step. It is on "running" for a while (like some seconds) and then it just stops and reset things like I quit kernel. I've experienced this kind of error? occasionally, but it was resolved just by rerunning, but this one never works. Can you also figure this out why? $\endgroup$
    – physgj
    Commented Nov 22, 2021 at 0:35
  • $\begingroup$ There is not enough info for me to make a recommendation. $\endgroup$
    – Bob Hanlon
    Commented Nov 22, 2021 at 0:58
  • $\begingroup$ I added my actual functions in the main post. If you can take a look at it. Thank you much anyway. $\endgroup$
    – physgj
    Commented Nov 22, 2021 at 1:33

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.