Message Boards Message Boards

How to improve the accuracy of FindRoot or NDSolveValue at extreme conditions in my equation system?

Mikoto Misaka
Posted 1 year ago

This is a question related to a DAE system originally posted on Mathematica Stack Exchange here. Since I haven't received an answer yet, I tried to come here for help. The equation system now I am eagering for an answer is

L = 10; l = L*(3/1000); DD = 0.01; 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}];

then 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. The shape of function fin[s] that I want

Actually for this equation system I have received two answers from to users in Mathematica Stack Exchange @AlexTrounev and @bbgodfrey. However, it is related to my eariler question and the parameters there are set to l = L*(3/10); DD = 1; . Now in this one I especially need the results in condition of smaller DD and l. I have found there codes work not very well at small DD and l. It seems there exists a bottleneck in function FindRoot and NDSolveValue at extreme conditions, for example, in this question they are set to

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

I’m curious about what causes such difficulty in solving these equations at this condition in function FindRoot and NDSolveValue. Are there any techniques that can help improve this situation? I am eagering for an answer.
POSTED BY: Mikoto Misaka
6 Replies
Sort By:
Mikoto Misaka
Mikoto Misaka
Posted 1 year ago

The codes I posted in this question may not be the best choices. I am open to any methods that can solve this problem!
POSTED BY: Mikoto Misaka

I tried the first few iterations of FindRoot, It brings about T=75, F=27300. The plot of m1 reveals that the function has a very high number of oscillations between 0 and L. I suppose this is difficult to manage for the numerical solver NDSolve. I am no expert, and I cannot say more.

The ":=" matter is about style of programming; it does not affect computing time.
POSTED BY: Gianluca Gorni
Mikoto Misaka
Mikoto Misaka
Posted 1 year ago

@Gianluca Gorni Thank you very much! The codes I posted in this question may not be the best choices. I am open to any methods that can solve this problem. For instance, I have also tried using the function NDSolveValue instead of ParametricNDSolve and FindRoot, here is the code

L = 10; l = L*(3/10); DD = 1; mu = 1;
sol = NDSolveValue[{\[Phi]1''[s] - (T[s] - F[s])*Sin[\[Phi]1[s]] == 
        0, \[Phi]2''[s] - T[s]*Sin[\[Phi]2[s]] == 0, T'[s] == 0, 
       F'[s] == 0, \[Phi]1[0] == 0, \[Phi]2'[L] == 
        0, \[Phi]1[l] == \[Phi]2[l] == 
        ArcTan[1/mu], \[Phi]2'[l] == \[Phi]1'[l], 
       x'[s] == Cos[\[Phi]1[s]], x[0] == 0, 
       x[l] == DD}, {\[Phi]1, \[Phi]2, T[L], F[L]}, {s, 0, L}, 
      Method -> {"Shooting", 
        "StartingInitialConditions" -> {\[Phi]2'[L] == 0}, 
        MaxIterations -> 1000}, AccuracyGoal -> 10, PrecisionGoal -> 10];fin[s_] := 
         Piecewise[{{sol[[1]][s], 0 <= s < l}, {sol[[2]][s], 
            l <= s <= L}}]; Plot[fin[s], {s, 0, L}]

The codes above work well when I choose l = L*(3/10); DD = 1;, but I actually need the answer for l = L*(3/1000); DD = 0.01;, and then they perform poorly. I am trying to adjust the StartingInitialConditions but I still cannot get a good result. I would appreciate any techniques that can improve this situation.
POSTED BY: Mikoto Misaka

What is your value of sol2? On my computer the FindRoot does not seem to converge.

A curiosity: why to you type

m1[s_, T_, F_] := 
         Evaluate[Evaluate[ϕ1[T, F] /. sol][s]];

with ":=" (delayed evaluation) followed by Evaluate, instead of

m1[s_, T_, F_] = \[Phi]1[T, F][s] /. sol;

with "=", immediate evaluation?
POSTED BY: Gianluca Gorni
Mikoto Misaka
Mikoto Misaka
Posted 1 year ago

@Gianluca Gorni Hi, I appreciate your reply: ). I am using this code just to illustrate my equation system. I also encountered difficulties in making the function FindRoot converge for my chosen parameters and initial seeds,

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

so I am seeking for any method that could make it converge at such chosen parameters. I would be grateful if you could share any techniques that could improve this situation. Moreover, I use ":=" to define m1 as a function, and to apply it in BCs. I am not sure if this is appropriate or not since I am new to Mathematica. Thanks again!
POSTED BY: Mikoto Misaka

POSTED BY: Gianluca Gorni
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract