3
$\begingroup$

Original:

I am running the following code to find the root of F , T and MM, here is an example

L = 10; l = L*(3/10); DD = 1; mu = 1; sol = 
 ParametricNDSolve[{\[Phi]1''[s] - (T - F)*Sin[\[Phi]1[s]] == 
    0, \[Phi]2''[s] - T*Sin[\[Phi]2[s]] == 0, \[Phi]1[0] == 
    0, \[Phi]2'[L] == MM, \[Phi]1[l] == \[Phi]2[l] == 
    ArcTan[1/mu]}, {\[Phi]1, \[Phi]2}, {s, 0, L}, {T, F, MM}]; 
m1[s_, T_, F_, MM_] := 
 Evaluate[Evaluate[\[Phi]1[T, F, MM] /. sol][s]]; 
m2[s_, T_, F_, MM_] := 
 Evaluate[Evaluate[\[Phi]2[T, F, MM] /. sol][s]]; 
dm1[s_, T_, F_, MM_] := Evaluate[D[m1[s, T, F, MM], s]]; 
dm2[s_, T_, F_, MM_] := Evaluate[D[m2[s, T, F, MM], s]]; 
BC1[T_?NumericQ, F_?NumericQ, MM_?NumericQ] := 
 NIntegrate[Cos[m1[s, T, F, MM]], {s, 0, l}] - DD; 
BC2[T_?NumericQ, F_?NumericQ, MM_?NumericQ] := -dm1[l, T, F, MM] - 
  T*(NIntegrate[Sin[m2[s, T, F, MM]], {s, l, L}]) + MM; 
BC3[T_?NumericQ, F_?NumericQ, MM_?NumericQ] := m2[L, T, F, MM]; sol2 =
  Monitor[
  FindRoot[{BC1[T, F, MM] == 0, BC2[T, F, MM] == 0, 
    BC3[T, F, MM] == 0}, {{T, 1}, {F, 1}, {MM, 1}}], {T, F, MM}]

I've noticed that it takes quite a long time to get a response, and there's no assurance of the solution's accuracy. It's worth noting that the parameter DD can range between 0.01 and 2, and this makes the code's performance decline further when DD is small.

I'm seeking advice on techniques that can improve the speed and accuracy of the FindRoot function. Thanks a lot!

Update1:

I have removed the variable MM because it seems unnecessary, and the boundary conditions have also been changed for simplicity.Here is my updated code:

  L = 10; l = L*(3/10); DD = 0.1; mu = 1;
         sol = 
             ParametricNDSolve[{ϕ1''[s] - (T - F)*Sin[ϕ1[s]] == 
                0, ϕ2''[s] - T*Sin[ϕ2[s]] == 0, ϕ1[0] == 
                0, ϕ2'[L] == 0, ϕ1[l] == ϕ2[l] == 
                ArcTan[1/mu]}, {ϕ1, ϕ2}, {s, 0, L}, {T, F}]; 
            m1[s_, T_, F_] := 
             Evaluate[Evaluate[ϕ1[T, F] /. sol][s]]; 
            m2[s_, T_, F_] := 
             Evaluate[Evaluate[ϕ2[T, F] /. sol][s]]; 
            dm1[s_, T_, F_] := Evaluate[D[m1[s, T, F], s]]; 
            dm2[s_, T_, F_] := Evaluate[D[m2[s, T, F], s]]; 
            BC1[T_?NumericQ, F_?NumericQ] := 
             NIntegrate[Cos[m1[s, T, F]], {s, 0, l}] - DD; 
            BC2[T_?NumericQ, F_?NumericQ] := -dm1[l, T, F] + dm2[l, T, F]; 
           sol2 =
              Monitor[
              FindRoot[{BC1[T, F] == 0, BC2[T, F] == 0}, {{T, 1}, {F, 1}}], {T, F}]

Moreover, define a Piecewise function which is the conbination of ϕ1 and ϕ2,

fin[s_] := Piecewise[{{ϕ1[s], 0 <= s < l}, {ϕ2[s], l <= s <= L}}];

Since for different initial values used in Findroot, we can get different ϕ1 and ϕ2. However, the shape of fin[s] shown in the following figure is actually my goal. enter image description here

Update2:

I especially need the results in condition of smaller DD and l. Such as

l = L*(2/100); DD = 0.1;

or even at a more extreme one

l = L*(2/1000); DD = 0.01;

I’m curious about what causes such difficulty in solving these equations at this condition in function FindRoot and NDSolve. Are there any techniques that can help improve this situation?

Improve this question
$\endgroup$
4
  • $\begingroup$ Are you sure that your problem has single solution? It seems that you don't use $\]Phi\1$` and `\[Phi]2` except for computing `\[Phi]1'` and `\[Phi]2'`. Why don't you instruct `ParametricNDSolve` to return `\[Phi]1'` and `\[Phi]2'` rather then `\[Phi]1` and `\[Phi]2` ? $\endgroup$
    – Igor Kotelnikov
    Commented Jul 8, 2023 at 14:00
  • $\begingroup$ @IgorKotelnikov. I appreciate your comment. I am new to Mathematica and I am not sure which method is the best. I would be grateful for any method that can solve this system of equations. I understand that the value of \[Phi]1 and \[Phi]2 is not unique, it will depend on the initial value of the solution and converge to different results. $\endgroup$
    – Mikoto
    Commented Jul 8, 2023 at 14:42
  • $\begingroup$ Hi. I am using similar shooting method to solve boundary-value problem. First of all, I would recommend explicitly restrict the range of expected values of free parameters in ParamentricNDSolve and instruct it to return only those functions which you need. Here is my case: ParametricNDSolveValue[ {ode$[k, z, p, {q, Mm, Rm}, \[CapitalLambda]] == 0, \[Phi][0] == 1, \[Phi]'[0] == 0} , {\[Phi]'[zR$], \[Phi][zR$]} , {z, 0, (1 + 0* 10^-7) zR$} , {{p, 0, pmx}, {\[CapitalLambda], 1, 500}} (*,Method->{"StiffnessSwitching"} ,Method->{"EquationSimplification"->"Residual"}*) ] $\endgroup$
    – Igor Kotelnikov
    Commented Jul 9, 2023 at 4:40
  • $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$
    – bbgodfrey
    Commented Jul 10, 2023 at 21:12

2 Answers 2

Reset to default
3
$\begingroup$

With g = T - F, the equations for ϕ1 and ϕ2 can be decoupled and solved sequentially. Consider ϕ1.

sL = 10; sm = sL*(3/10); dd = 1; mu = 1;

and

{ϕ1''[s] - g*Sin[ϕ1[s]] == 0, ϕ1[0] == 0, ϕ1[sm] == ArcTan[1/mu]};

(I have changed the names of some parameters to avoid using capital letters and for clarity.)

There remains the BC1 boundary condition, which requires that the integral of Cos[ϕ1[s]]over {s, 0, sm} be equal to dd. Introduce the new variable d'[s] = Cos[ϕ1[s]]. Then this boundary condition becomes simply d[sm] - d[0] == dd.

Finally, determining the eigenvalue g can be accomplished using the procedure described here. Combining all these yields

sol1 = NDSolveValue[{ϕ1''[s] - g[s]*Sin[ϕ1[s]] == 0, g'[s] == 0, 
    d'[s] == Cos[ϕ1[s]], ϕ1[0] == 0, ϕ1[sm] == ArcTan[1/mu], 
    d[0] == 0, d[sm] == dd}, {ϕ1[s], ϕ1'[s], g[sm]}, {s, 0, sL}];

sol1[[3]
(* -1.13517 *)
NIntegrate[Cos[sol1[[1]]], {s, 0, sm}]
(* 1. *)

where the first value is g and the second verifies that the integral of Cos[ϕ1[s]] is indeed dd. And,

Plot[Evaluate@sol1[[;; 2]], {s, 0, sL}, AxesLabel -> {s, "ϕ1,ϕ1'"},
    PlotLabels -> Placed[{"ϕ1", "ϕ1'"}, {{Scaled[1], Before}}], 
    LabelStyle -> {12, Bold, Black}]

enter image description here

Next, consider ϕ2.

{ϕ2''[s] - t*Sin[ϕ2[s]] == 0, ϕ2[sm] == ArcTan[1/mu, ϕ2'[sL] == mm];

together with BC2 and BC3. The integral in BC2 is solved analytically by integrating the ϕ2 ODE once to obtain ϕ2[sL] - ϕ2[sm], obviating the need to introduce an auxiliary variable similar d for ϕ1. The two remaining terms in this boundary condition are given by ϕ2'[sL] (mm from the third equation above) and ϕ1'[sm]. Combined, ϕ2'[sm] == (sol1[[2]] /. s -> sm). In other words, ϕ1 and ϕ2 are tangent at sm. Finally, BC3 is equivalent to ϕ2[sL] == 0. In total,

sol2 = NDSolveValue[{ϕ2''[s] - t[s]*Sin[ϕ2[s]] == 0, t'[s] == 0, 
    ϕ2'[sm] == (sol1[[2]] /. s -> sm), ϕ2[sm] == ArcTan[1/mu],
    ϕ2[sL] == 0}, {ϕ2[s], ϕ2'[s], t[sm]}, {s, 0, sL}, 
    InitialSeeding -> t[sL] == -1/2];

sol2[[3]]
(* -0.597091 *)

is the value of the eigenvalue t here. Note that many values of t satisfy the ϕ2 system. I chose InitialSeeding -> t[sL] == -1 to bias the calculation toward a solution resembling ϕ1 above.

Plot[Evaluate@sol2[[;; 2]], {s, 0, sL}, AxesLabel -> {s, "ϕ2,ϕ2'"},
    PlotLabels -> Placed[{"ϕ2", "ϕ2'"}, {{Scaled[1], Before}}],
    LabelStyle -> {12, Bold, Black}]

enter image description here

Overlay ϕ1 and ϕ2 for completeness. As expected, the curves are tangent at sm.

Plot[{sol1[[1]], sol2[[1]]}, {s, 0, sL}, AxesLabel -> {s, "ϕ1,ϕ2"},
    PlotLabels -> Placed[{"ϕ1", "ϕ2"}, {{Scaled[1], Before}}],
    LabelStyle -> {12, Bold, Black}]

enter image description here

Incidentally, the solution obtained by Alex Trounev can be obtained by using InitialSeeding -> t[sL] == -1/2.

Improve this answer
$\endgroup$
6
  • $\begingroup$ @ bbgodfrey Thank you very much for your help! However, I encountered some difficulties with your method when the parameter DD is small. For example, when DD=0.1, Mathematica does not find the accurate solution of T and ϕ2. Do you have any suggestions on how to overcome this inaccuracy? By the way, I have eliminated the variable MM to make the question clearer. $\endgroup$
    – Mikoto
    Commented Jul 16, 2023 at 16:01
  • $\begingroup$ With my code above, I obtained a result for dd =1/10 without difficulty. It looks like the result in my answer, but with ϕ1 and ϕ2 overlapping more closely. I see that you have added a curve to your question. How did you obtain this result, and what are the associated values of your input parameters and of F and T?. $\endgroup$
    – bbgodfrey
    Commented Jul 17, 2023 at 4:30
  • $\begingroup$ Thanks again! To get the curve I added to the question, I choose the initial parameters：sL = 10; sm = sL*(2/10); dd = 1; mu = 1. For sol1, I have added InitialSeeding -> {[Phi]1'[0] == 2}, thus getting g[sm]=-2.08538. For sol2, here is my following code: sol2 = NDSolveValue[{[Phi]2''[s] - t[s]*Sin[[Phi]2[s]] == 0, t'[s] == 0, [Phi]2'[sm] == (Evaluate[sol1[[2]]][s] /. s -> sm), [Phi]2[sm] == ArcTan[1/mu], [Phi]2'[sL] == 0}, {[Phi]2, [Phi]2', t[sm]}, {s, sm, sL}, InitialSeeding -> {[Phi]2[sL] == 0}]. The boundary conditions have been changed since I have deleted the variable MM. $\endgroup$
    – Mikoto
    Commented Jul 17, 2023 at 6:13
  • $\begingroup$ You also seem to have replaced ϕ2[sL] == 0 (the equivalent of BC3 == 0) by ϕ2'[sL] == 0. $\endgroup$
    – bbgodfrey
    Commented Jul 17, 2023 at 16:44
  • $\begingroup$ I agree with the comment below by @AlexTrounev that you should not change a question that already has received answers. If you have a new question, post it as such. $\endgroup$
    – bbgodfrey
    Commented Jul 18, 2023 at 3:38
3
$\begingroup$

We can add BC1,BC2,BC3 to the system of equations and solve it as DAE system. Unfortunately NDSolve can't solve this type DAE. Therefore we use our favorable the Euler wavelets collocation method described here, here, and here. For this we map system on $0\le s\le1$ using scale L, we have

L = 10; l = (3/
   10); DD = 1; mu = 1; {ϕ1''[s] - (T - F)*Sin[ϕ1[s]] == 
  0, ϕ2''[s] - T*Sin[ϕ2[s]] == 0, f1'[s] == Cos[ϕ1[s]],
  f2'[s] == Sin[ϕ2[s]], ϕ1[0] == 0, ϕ2'[L] == 
  MM, ϕ1[l] == ArcTan[1/mu], ϕ2[l] == ArcTan[1/mu], 
 f1[0] == 0, 
 f2[0] == 0}; eqs = {p1''[s]/L^2 - (T - F)*Sin[p1[s]] == 0, 
  p2''[s]/L^2 - T*Sin[p2[s]] == 0, f1'[s]/L == Cos[p1[s]], 
  f2'[s]/L == Sin[p2[s]]}; bc = {p1[0] == 0, p2'[1]/L == MM, 
  p1[l] == ArcTan[1/mu], p2[l] == ArcTan[1/mu], f1[0] == 0, 
  f2[0] == 0, p2[1] == 0, MM - p1'[l]/L - T (f2[1] - f2[l]) == 0, 
  f1[l] - DD == 0}; rule = {p1'' -> ddp1, p1' -> dp1, p2'' -> ddp2, 
  p2' -> dp2, f1' -> df1, f2' -> df2};
eqs1 = eqs /. rule; bc1 = bc /. rule;
UE[m_, t_] := EulerE[m, t];
psi[k_, n_, m_, t_] := 
  Piecewise[{{2^(k/2) UE[m, 2^k t - 2 n + 1], (n - 1)/2^(k - 1) <= t <
       n/2^(k - 1)}, {0, True}}];
PsiE[k_, M_, t_] := 
 Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 6; M0 = 7; With[{k = k0, M = M0}, 
 nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; ycol = 
 Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk = 
 With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 = 
 With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y;

P1 = Array[pp1, {nn}]; P2 = Array[pp2, {nn}]; F1 = 
 Array[ff1, {nn}]; F2 = Array[ff2, {nn}];
ddp1[s_] := P1 . Psi[s]; dp1[s_] := P1 . int1[s] + a0; 
p1[s_] := P1 . int2[s] + a0 s + a1;
ddp2[s_] := P2 . Psi[s]; dp2[s_] := P2 . int1[s] + b0; 
p2[s_] := P2 . int2[s] + b0 s + b1;
df1[s_] := F1 . Psi[s]; f1[s_] := F1 . int1[s] + c0; 
df2[s_] := F2 . Psi[s]; f2[s_] := F2 . int1[s] + d0;
eqn = Table[eqs1, {s, ycol}];
eqsAll = Join[Flatten[eqn], bc];
var = Join[{T, F, MM, a0, a1, b0, b1, c0, d0}, P1, P2, F1, F2];

Numerical solution (note that we use FindRoot one time only!)

sol1 = FindRoot[eqsAll, Table[{var[[i]], 1/10}, {i, Length[var]}]];

Visualization

plot1= Plot[Evaluate[{p1[s/L], p2[s/L]} /. sol1], {s, 0, L}, 
     PlotLegends -> {ϕ1, ϕ2}]

Figure1

Now we can solve this problem with NDSolve using parameters

rule1 = Take[sol1, 3]

(*Out[]= {T -> -0.597091, F -> 0.53808, MM -> 1.44889}*)
L = 10; l = 
 L*(3/10); DD = 1; mu = 1; 
solq = NDSolveValue[{ϕ1''[s] - (T - F)*Sin[ϕ1[s]] == 
     0, ϕ2''[s] - T*Sin[ϕ2[s]] == 0, 
    f1'[s] == Cos[ϕ1[s]], 
    f2'[s] == Sin[ϕ2[s]], ϕ1[0] == 0, ϕ2'[L] == 
     MM, ϕ2[L] == 0, f1[0] == 0, f2[0] == 0, f1[l] == DD} /. 
   rule1, {ϕ1, ϕ2, ϕ1', ϕ2', f1, f2}, {s, 0, L}, 
  Method -> {"Shooting", 
    "ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"}}];

Visualization together with plot1

Show[plot1, Plot[Table[solq[[i]][s], {i, 2}] // Evaluate, {s, 0, L}, 
  PlotLegends -> {ϕ1, ϕ2}, 
  PlotStyle -> {{Red, Dashed}, {Blue, Dashed}}]]

Figure 2

As we can see two solutions consider well. Finally note, that we update boundary conditions using BC1,BC3 instead of ϕ1[l] == ϕ2[l] == ArcTan[1/mu] since NDSolve can't handle boundary conditions in the intermediate point s=l with high accuracy. We can check that solq satisfiers these equations with error of

{solq[[1]][l] - ArcTan[1/mu], solq[[2]][l] - ArcTan[1/mu]}

(*Out[]= {1.06451*10^-8, -2.84661*10^-7}*)

Update 1. Solution for updated problem

L = 10; l = (3/10); DD = 1; mu = 1; {ϕ1''[s] - (T - F)*Sin[ϕ1[s]] == 
  0, ϕ2''[s] - T*Sin[ϕ2[s]] == 0, f1'[s] == Cos[ϕ1[s]],
  f2'[s] == Sin[ϕ2[s]], ϕ1[0] == 0, ϕ2'[L] == 
  MM, ϕ1[l] == ArcTan[1/mu], ϕ2[l] == ArcTan[1/mu], 
 f1[0] == 0, 
 f2[0] == 0}; eqs = {p1''[s]/L^2 - (T - F)*Sin[p1[s]] == 0, 
  p2''[s]/L^2 - T*Sin[p2[s]] == 0, 
  f1'[s]/L == Cos[p1[s]]}; bc = {p1[0] == 0, p1[l] == ArcTan[1/mu], 
  p2[l] == ArcTan[1/mu], f1[0] == 0, p2[1] == 0, 
  p2'[l]/L - p1'[l]/L == 0, f1[l] - DD == 0}; rule = {p1'' -> ddp1, 
  p1' -> dp1, p2'' -> ddp2, p2' -> dp2, f1' -> df1};
eqs1 = eqs /. rule; bc1 = bc /. rule;
UE[m_, t_] := EulerE[m, t];
psi[k_, n_, m_, t_] := 
  Piecewise[{{2^(k/2) UE[m, 2^k t - 2 n + 1], (n - 1)/2^(k - 1) <= t <
       n/2^(k - 1)}, {0, True}}];
PsiE[k_, M_, t_] := 
 Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 6; M0 = 7; With[{k = k0, M = M0}, 
 nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; ycol = 
 Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk = 
 With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 = 
 With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y;

P1 = Array[pp1, {nn}]; P2 = Array[pp2, {nn}]; F1 = Array[ff1, {nn}];
ddp1[s_] := P1 . Psi[s]; dp1[s_] := P1 . int1[s] + a0;
p1[s_] := P1 . int2[s] + a0 s + a1;
ddp2[s_] := P2 . Psi[s]; dp2[s_] := P2 . int1[s] + b0;
p2[s_] := P2 . int2[s] + b0 s + b1;
df1[s_] := F1 . Psi[s]; f1[s_] := F1 . int1[s] + c0;
eqn = Table[eqs1, {s, ycol}];
eqsAll = Join[Flatten[eqn], bc];
sol2 = FindRoot[eqsAll, Table[{var[[i]], 1/10}, {i, Length[var]}]];

plot2 = Plot[Evaluate[{p1[s/L], p2[s/L]} /. sol2], {s, 0, L}, 
  PlotLegends -> {ϕ1, ϕ2}, PlotRange -> All]

fin[s_] := 
  Piecewise[{{p1[s/L] /. sol2, 0 <= s < l L}, {p2[s/L] /. sol2, 
     l L <= s <= L}}];

Plot[Evaluate[fin[s]], {s, 0, L}]

Figure 3

Update 2 In the code above we put l = (2/10); DD = 1/10;, results shown below Figure 4

Improve this answer
$\endgroup$
11
  • $\begingroup$ @ Alex Trounev. Thanks a lot for your help! I believe your method could solve my problem completely. However, I'm still learning your code since I am new to the Euler wavelets collocation method, and I would like to know how to adjust the initial seeds to get different ϕ1 and ϕ2 through your method. By the way, I have removed the variable MM to simplify the question. $\endgroup$
    – Mikoto
    Commented Jul 16, 2023 at 16:08
  • 1
    $\begingroup$ @Mikoto Did you offer a new problem instead of old one? It seems to be against rule of this forum as I understand. Normally we use update for the old question with new thoughts. $\endgroup$
    – Alex Trounev
    Commented Jul 17, 2023 at 7:20
  • $\begingroup$ @ Alex Trounev.I tried to make the problem clearer by simplifying it instead of offering a new one. I removed the variable MM because it seemed unnecessary and it made the problem more difficult. I also shared some of the results I am aiming for, but I don’t know how to pick the initial seeds to get them. I don't know whether it will against the rules this forum since I am new to here. Please forgive me if this modification makes you confused. $\endgroup$
    – Mikoto
    Commented Jul 17, 2023 at 7:44
  • 1
    $\begingroup$ @Mikoto It is not confusing me, it is not clear for readers who will try to understand to which task I gave my solution. :) $\endgroup$
    – Alex Trounev
    Commented Jul 17, 2023 at 9:57
  • $\begingroup$ I apologize for my mistake. Now I have restored the question and posted an updated one, which won't confuse any readers. $\endgroup$
    – Mikoto
    Commented Jul 18, 2023 at 6:37

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.