1
$\begingroup$

I am obtaining two different answers for a curve when solving the same exact system of differential equations when using NDSolve vs when using ParametricNDSolveValue. Here it is when using NDSolve:

h[r_] = 10^(-2); 
boundary1 = 0.5; 
boundary2 = 1; 
\[Mu]1 = 8*10^4; 
\[Mu]2 = 8*10^4; 
p = 1300; 
w = 1.48323; 
nsolution1 = NDSolve[{lt[r] + r*Derivative[1][lt][r] == lr[r]*Cos[\[Beta][r]], p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]1*h[r]*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h[r]*r*Derivative[1][\[Beta]][r]*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - 
        \[Mu]1*h[r]*r*Sin[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) == 0, (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]1*h[r]*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + 
        \[Mu]1*h[r]*r*Derivative[1][\[Beta]][r]*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h[r]*r*Cos[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]1*h[r]*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0, lt[0.001] == w, 
      lr[0.001] == w, \[Beta][0.001] == 0.001}, {lr, lt, \[Beta]}, {r, 0.001, boundary1}][[1]]; 
radialstretch[(r_)?NumericQ] := lr[r] /. nsolution1
thetastretch[(r_)?NumericQ] := lt[r] /. nsolution1
angle[(r_)?NumericQ] := \[Beta][r] /. nsolution1
radialstretchderivative[(r_)?NumericQ] := D[lr[u] /. nsolution1, u] /. u -> r
thetastretchderivative[(r_)?NumericQ] := D[lt[u] /. nsolution1, u] /. u -> r
R[(r_)?NumericQ] := r*thetastretch[r]
Rderivative[(r_)?NumericQ] := thetastretch[r] + r*thetastretchderivative[r]
f1[(r_)?NumericQ] := -(radialstretch[r]^2 - Rderivative[r]^2)^(1/2)
Z[(r_)?NumericQ] := NIntegrate[Re[f1[t]], {t, 0.001, r}]
p1 := ParametricPlot[{R[r], Z[r] + 1}, {r, 0.001, boundary1}, PlotStyle -> Green]
g[r_] := (1/Sqrt[2])*((\[Mu]1/\[Mu]2)*(radialstretch[r]^2 - 1/(radialstretch[r]^2*thetastretch[r]^2)) + ((\[Mu]1 - radialstretch[r]^4*thetastretch[r]^2*\[Mu]1)^2 + 4*radialstretch[r]^4*thetastretch[r]^2*\[Mu]2^2)^(1/2)/(radialstretch[r]^2*thetastretch[r]^2*\[Mu]2))^(1/2)
nsolution2 = Quiet[NDSolve[{lt[r] + r*Derivative[1][lt][r] == lr[r]*Cos[\[Beta][r]], p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]2*h[r]*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h[r]*r*Derivative[1][\[Beta]][r]*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - 
         \[Mu]2*h[r]*r*Sin[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) == 0, (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]2*h[r]*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + 
         \[Mu]2*h[r]*r*Derivative[1][\[Beta]][r]*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h[r]*r*Cos[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]2*h[r]*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0, 
       lt[boundary1] == thetastretch[boundary1], lr[boundary1] == g[boundary1], \[Beta][boundary1] == angle[boundary1]}, {lr, lt, \[Beta]}, {r, boundary1, boundary2}][[1]]]; 
radialstretch2[(r_)?NumericQ] := Quiet[lr[r] /. nsolution2]
thetastretch2[(r_)?NumericQ] := Quiet[lt[r] /. nsolution2]
angle2[(r_)?NumericQ] := Quiet[\[Beta][r] /. nsolution2]
radialstretchderivative2[(r_)?NumericQ] := Quiet[D[lr[u] /. nsolution2, u] /. u -> r]
thetastretchderivative2[(r_)?NumericQ] := Quiet[D[lt[u] /. nsolution2, u] /. u -> r]
R2[(r_)?NumericQ] := r*thetastretch2[r]
Rderivative2[(r_)?NumericQ] := thetastretch2[r] + r*thetastretchderivative2[r]
f2[(r_)?NumericQ] := Quiet[-(radialstretch2[r]^2 - Rderivative2[r]^2)^(1/2)]
Z2[(r_)?NumericQ] := Quiet[NIntegrate[Re[f2[t]], {t, boundary1, r}]]
sol2 = Solve[Z2[boundary1] + x2 == Z[boundary1], Plus[x2]]; 
x2 = x2 /. sol2[[1]]; 
y2 := 1 + x2; 
p2 := ParametricPlot[{R2[r], Z2[r] + y2}, {r, boundary1, boundary2}, PlotStyle -> Red]
Show[p1, p2, PlotRange -> All]

curve shape when using NDSolve

Now, here it is solving it using ParametricNDSolveValue:

h = 1/10^2; 
\[Mu]1 = 8*10^4; 
\[Mu]2 = 8*10^4; 
p = 1300; 
boundary1 = 0.5; 
boundary2 = 1; 
nsolution1 = ParametricNDSolve[{lt[r] + r*Derivative[1][lt][r] == lr[r]*Cos[\[Beta][r]], p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]1*h*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h*r*Derivative[1][\[Beta]][r]*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - 
       \[Mu]1*h*r*Sin[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) == 0, 
     (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]1*h*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + \[Mu]1*h*r*Derivative[1][\[Beta]][r]*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h*r*Cos[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + 
         (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]1*h*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0, lt[0.001] == w, lr[0.001] == w, \[Beta][0.001] == 0.001}, {lr, lt, \[Beta]}, {r, 0.001, boundary1}, {w}, AccuracyGoal -> 10, PrecisionGoal -> 10]; 
radialstretch1[(w_)?NumericQ, (r_)?NumericQ] := lr[w][r] /. nsolution1
thetastretch1[(w_)?NumericQ, (r_)?NumericQ] := lt[w][r] /. nsolution1
angle1[(w_)?NumericQ, (r_)?NumericQ] := \[Beta][w][r] /. nsolution1
radialstretchderivative1[(w_)?NumericQ, (r_)?NumericQ] := D[lr[w][t] /. nsolution1, t] /. t -> r
thetastretchderivative1[(w_)?NumericQ, (r_)?NumericQ] := D[lt[w][t] /. nsolution1, t] /. t -> r
R1[(w_)?NumericQ, (r_)?NumericQ] := r*thetastretch1[w, r]
Rderivative1[(w_)?NumericQ, (r_)?NumericQ] := thetastretch1[w, r] + r*thetastretchderivative1[w, r]
f1[(w_)?NumericQ, (r_)?NumericQ] := -Sqrt[radialstretch1[w, r]^2 - Rderivative1[w, r]^2]
Z1[(w_)?NumericQ, (r_)?NumericQ] := NIntegrate[f1[w, r], {t, 0.001, r}]
g1[(w_)?NumericQ, (r_)?NumericQ] := (1/Sqrt[2])*((\[Mu]1/\[Mu]2)*(radialstretch1[w, r]^2 - 1/(radialstretch1[w, r]^2*thetastretch1[w, r]^2)) + ((\[Mu]1 - radialstretch1[w, r]^4*thetastretch1[w, r]^2*\[Mu]1)^2 + 4*radialstretch1[w, r]^4*thetastretch1[w, r]^2*\[Mu]2^2)^(1/2)/
      (radialstretch1[w, r]^2*thetastretch1[w, r]^2*\[Mu]2))^(1/2)
nsolution2 = ParametricNDSolve[{lt[r] + r*Derivative[1][lt][r] == lr[r]*Cos[\[Beta][r]], p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]2*h*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h*r*Derivative[1][\[Beta]][r]*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - 
       \[Mu]2*h*r*Sin[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) == 0, 
     (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]2*h*Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + \[Mu]2*h*r*Derivative[1][\[Beta]][r]*Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h*r*Cos[\[Beta][r]]*(Derivative[1][lr][r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + 
         (2*Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]2*h*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0, lt[boundary1] == thetastretch1[w, boundary1], lr[boundary1] == g1[w, boundary1], \[Beta][boundary1] == angle1[w, boundary1]}, {lr, lt, \[Beta]}, {r, boundary1, boundary2}, {w}, AccuracyGoal -> 10, 
    PrecisionGoal -> 10]; 
radialstretch2[(w_)?NumericQ, (r_)?NumericQ] := lr[w][r] /. nsolution2
thetastretch2[(w_)?NumericQ, (r_)?NumericQ] := lt[w][r] /. nsolution2
angle2[(w_)?NumericQ, (r_)?NumericQ] := \[Beta][w][r] /. nsolution2
radialstretchderivative2[(w_)?NumericQ, (r_)?NumericQ] := D[lr[w][t] /. nsolution2, t] /. t -> r
thetastretchderivative2[(w_)?NumericQ, (r_)?NumericQ] := D[lt[w][t] /. nsolution2, t] /. t -> r
R2[(w_)?NumericQ, (r_)?NumericQ] := r*thetastretch2[w, r]
Rderivative2[(w_)?NumericQ, (r_)?NumericQ] := thetastretch2[w, r] + r*thetastretchderivative2[w, r]
f2[(w_)?NumericQ, (r_)?NumericQ] := -Sqrt[radialstretch2[w, r]^2 - Rderivative2[w, r]^2]
Z2[(w_)?NumericQ, (r_)?NumericQ] := NIntegrate[f2[w, r], {t, boundary1, r}]
p1 := ParametricPlot[{R1[1.48323, r], Z1[1.48323, r] + 1}, {r, 0.001, boundary1}, PlotStyle -> Green]
Clear[x2]; 
Clear[y2]; 
sol2 = Solve[Z2[1.48323, boundary1] + x2 == Z1[1.48323, boundary1], Plus[x2]]; 
x2 = x2 /. sol2[[1]]; 
y2 = 1 + x2; 
p2 := ParametricPlot[{R2[1.48323, r], Z2[1.48323, r] + y2}, {r, boundary1, boundary2}, PlotStyle -> Red]
Show[p1, p2, PlotRange -> All]

curve shape when using ParametricNDSolveValue

As you can see, the equations are exactly the same, the parameters used are exactly the same, however, the results are different. I would love to use ParametricNDSolveValue because it is much more convenient for my code, which is much longer than this, and it works much faster, but I need to first know why it is giving me a different answer. If anyone has any insight, please let me know.

Improve this question
$\endgroup$
2
  • $\begingroup$ There are many typos in the code, for example, h[r] is not defined for nsolution1. $\endgroup$
    – Alex Trounev
    Commented Jan 29, 2023 at 5:59
  • $\begingroup$ It was in my notebook, I somehow missed copying it here. I edited it now. $\endgroup$
    – juv95
    Commented Jan 29, 2023 at 7:03

1 Answer 1

Reset to default
1
$\begingroup$

Solutions computed with NDSolve and ParametricNDSolve are the same as it follows from the test

Clear["`*"]

boundary1 = 0.5;
boundary2 = 1;
\[Mu]1 = 8*10^4;
\[Mu]2 = 8*10^4;
p = 1300;
h = 1/10^2;

nsolutionp1 = 
 ParametricNDSolve[{lt[r] + r*Derivative[1][lt][r] == 
    lr[r]*Cos[\[Beta][r]], 
   p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]1*h*
      Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h*r*
      Derivative[1][\[Beta]][r]*
      Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h*r*
      Sin[\[Beta][
        r]]*(Derivative[1][lr][
         r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
           Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) == 
    0, (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]1*h*
      Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + \[Mu]1*h*r*
      Derivative[1][\[Beta]][r]*
      Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h*r*
      Cos[\[Beta][
        r]]*(Derivative[1][lr][
         r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
           Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]1*
      h*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0, lt[0.001] == w, 
   lr[0.001] == w, \[Beta][0.001] == 0.001}, {lr, lt, \[Beta]}, {r, 
   0.001, boundary1}, {w}]; 
radialstretch1[(w_)?NumericQ, (r_)?NumericQ] := 
 lr[w][r] /. nsolutionp1;
thetastretch1[(w_)?NumericQ, (r_)?NumericQ] := lt[w][r] /. nsolutionp1;
angle1[(w_)?NumericQ, (r_)?NumericQ] := \[Beta][w][r] /. nsolutionp1;
radialstretchderivative1[(w_)?NumericQ, (r_)?NumericQ] := 
  D[lr[w][t], t] /. t -> r /. nsolutionp1;
thetastretchderivative1[(w_)?NumericQ, (r_)?NumericQ] := 
  D[lt[w][t], t] /. t -> r /. nsolutionp1;
R1[(w_)?NumericQ, (r_)?NumericQ] := r*thetastretch1[w, r];
Rderivative1[(w_)?NumericQ, (r_)?NumericQ] := 
  thetastretch1[w, r] + r*thetastretchderivative1[w, r];
f1[(w_)?NumericQ, (r_)?NumericQ] := -Sqrt[
    radialstretch1[w, r]^2 - Rderivative1[w, r]^2];
Z1[(w_)?NumericQ, (r_)?NumericQ] := 
  NIntegrate[f1[w, r], {t, 0.001, r}];
g1[(w_)?NumericQ, (r_)?
   NumericQ] := (1/
    Sqrt[2])*((\[Mu]1/\[Mu]2)*(radialstretch1[w, r]^2 - 
        1/(radialstretch1[w, r]^2*
           thetastretch1[w, r]^2)) + ((\[Mu]1 - 
            radialstretch1[w, r]^4*thetastretch1[w, r]^2*\[Mu]1)^2 + 
         4*radialstretch1[w, r]^4*thetastretch1[w, r]^2*\[Mu]2^2)^(1/
         2)/(radialstretch1[w, r]^2*thetastretch1[w, r]^2*\[Mu]2))^(1/
     2); nsolutionp2 = 
 ParametricNDSolve[{lt[r] + r*Derivative[1][lt][r] == 
    lr[r]*Cos[\[Beta][r]], 
   p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]2*h*
      Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h*r*
      Derivative[1][\[Beta]][r]*
      Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h*r*
      Sin[\[Beta][
        r]]*(Derivative[1][lr][
         r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
           Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) == 
    0, (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]2*h*
      Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + \[Mu]2*h*r*
      Derivative[1][\[Beta]][r]*
      Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h*r*
      Cos[\[Beta][
        r]]*(Derivative[1][lr][
         r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
           Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]2*
      h*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0, 
   lt[boundary1] == thetastretch1[w, boundary1], 
   lr[boundary1] == g1[w, boundary1], \[Beta][boundary1] == 
    angle1[w, boundary1]}, {lr, lt, \[Beta]}, {r, boundary1, 
   boundary2}, {w}];

To test solutions we plot

{test1 = 
  Plot[Evaluate[{lr[w][r], lt[w][r], \[Beta][w][r]} /. 
      w -> 1.48323 /. nsolutionp1], {r, 0.001, boundary1}, 
   PlotLegends -> {lr, lt, \[Beta]}], 
 test2 = Plot[
   Evaluate[{lr[w][r], lt[w][r], \[Beta][w][r]} /. w -> 1.48323 /. 
     nsolutionp2], {r, boundary1, boundary2}, 
   PlotLegends -> {lr, lt, \[Beta]}]}

Figure 1

NDSolve code

Clear[{h, w, f1}]

boundary1 = 0.5;
boundary2 = 1;
\[Mu]1 = 8*10^4;
\[Mu]2 = 8*10^4;
p = 1300;
w = 1.48323; h[r_] := 1/10^2;

nsolution1 = 
  NDSolve[{lt[r] + r*Derivative[1][lt][r] == lr[r]*Cos[\[Beta][r]], 
     p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]1*h[r]*
        Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h[r]*r*
        Derivative[1][\[Beta]][r]*
        Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h[r]*r*
        Sin[\[Beta][
          r]]*(Derivative[1][lr][
           r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
             Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) == 
      0, (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]1*h[r]*
        Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + \[Mu]1*h[r]*r*
        Derivative[1][\[Beta]][r]*
        Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]1*h[r]*r*
        Cos[\[Beta][
          r]]*(Derivative[1][lr][
           r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
             Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]1*
        h[r]*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0, lt[0.001] == w, 
     lr[0.001] == w, \[Beta][0.001] == 0.001}, {lr, lt, \[Beta]}, {r, 
     0.001, boundary1}][[1]];
radialstretch[(r_)?NumericQ] := lr[r] /. nsolution1
thetastretch[(r_)?NumericQ] := lt[r] /. nsolution1
angle[(r_)?NumericQ] := \[Beta][r] /. nsolution1
radialstretchderivative[(r_)?NumericQ] := 
 D[lr[u] /. nsolution1, u] /. u -> r
thetastretchderivative[(r_)?NumericQ] := 
 D[lt[u] /. nsolution1, u] /. u -> r
R[(r_)?NumericQ] := r*thetastretch[r]
Rderivative[(r_)?NumericQ] := 
 thetastretch[r] + r*thetastretchderivative[r]
f1[(r_)?NumericQ] := -(radialstretch[r]^2 - Rderivative[r]^2)^(1/2)
Z[(r_)?NumericQ] := NIntegrate[Re[f1[t]], {t, 0.001, r}]
p1 := ParametricPlot[{R[r], Z[r] + 1}, {r, 0.001, boundary1}, 
  PlotStyle -> Green]
g[r_] := (1/
    Sqrt[2])*((\[Mu]1/\[Mu]2)*(radialstretch[r]^2 - 
        1/(radialstretch[r]^2*
           thetastretch[r]^2)) + ((\[Mu]1 - 
            radialstretch[r]^4*thetastretch[r]^2*\[Mu]1)^2 + 
         4*radialstretch[r]^4*thetastretch[r]^2*\[Mu]2^2)^(1/
         2)/(radialstretch[r]^2*thetastretch[r]^2*\[Mu]2))^(1/2)
nsolution2 = 
  NDSolve[{lt[r] + r*Derivative[1][lt][r] == lr[r]*Cos[\[Beta][r]], 
     p*r*lt[r]*lr[r]*Cos[\[Beta][r]] - \[Mu]2*h[r]*
        Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h[r]*r*
        Derivative[1][\[Beta]][r]*
        Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h[r]*r*
        Sin[\[Beta][
          r]]*(Derivative[1][lr][
           r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
             Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) == 
      0, (-p)*r*lr[r]*lt[r]*Sin[\[Beta][r]] - \[Mu]2*h[r]*
        Cos[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) + \[Mu]2*h[r]*r*
        Derivative[1][\[Beta]][r]*
        Sin[\[Beta][r]]*(lr[r] - 1/(lr[r]^3*lt[r]^2)) - \[Mu]2*h[r]*r*
        Cos[\[Beta][
          r]]*(Derivative[1][lr][
           r] + (3*Derivative[1][lr][r])/(lr[r]^4*lt[r]^2) + (2*
             Derivative[1][lt][r])/(lr[r]^3*lt[r]^3)) + \[Mu]2*
        h[r]*(lt[r] - 1/(lr[r]^2*lt[r]^3)) == 0, 
     lt[boundary1] == thetastretch[boundary1], 
     lr[boundary1] == g[boundary1], \[Beta][boundary1] == 
      angle[boundary1]}, {lr, lt, \[Beta]}, {r, boundary1, 
     boundary2}][[1]];
radialstretch2[(r_)?NumericQ] := Quiet[lr[r] /. nsolution2]
thetastretch2[(r_)?NumericQ] := Quiet[lt[r] /. nsolution2]
angle2[(r_)?NumericQ] := Quiet[\[Beta][r] /. nsolution2]
radialstretchderivative2[(r_)?NumericQ] := 
 Quiet[D[lr[u] /. nsolution2, u] /. u -> r]
thetastretchderivative2[(r_)?NumericQ] := 
 Quiet[D[lt[u] /. nsolution2, u] /. u -> r]
R2[(r_)?NumericQ] := r*thetastretch2[r]
Rderivative2[(r_)?NumericQ] := 
 thetastretch2[r] + r*thetastretchderivative2[r]
f2[(r_)?NumericQ] := 
 Quiet[-(radialstretch2[r]^2 - Rderivative2[r]^2)^(1/2)]
Z2[(r_)?NumericQ] := Quiet[NIntegrate[Re[f2[t]], {t, boundary1, r}]]
sol2 = Solve[Z2[boundary1] + x2 == Z[boundary1], x2];
x2 = x2 /. sol2[[1]];
y2 := 1 + x2;
p2 := ParametricPlot[{R2[r], Z2[r] + y2}, {r, boundary1, boundary2}, 
  PlotStyle -> Red]

Now we can compare two solutions as follows

{Show[test1, 
  Plot[Evaluate[{lr[r], lt[r], \[Beta][r]} /. nsolution1], {r, 0.001, 
    boundary1}, PlotLegends -> {lr, lt, \[Beta]}, 
   PlotStyle -> {{Black, Dashed}, {Black, Dashed}, {Black, 
      Dashed}}]],
 Show[test2, 
  Plot[Evaluate[{lr[r], lt[r], \[Beta][r]} /. nsolution2], {r, 
    boundary1, boundary2}, PlotLegends -> {lr, lt, \[Beta]}, 
   PlotStyle -> {{Black, Dashed}, {Black, Dashed}, {Black, Dashed}}]]}

Figure 2

Therefore, all discrepancies coming from definitions p1,p2 and other functions. For example, Z and Z1 are differ, while definitions f1 are not differ
Figure 3

We should look what NIntegrate doing in two cases. Finally we found out typos in definitions Z1[w, r], Z2[w,r]

Z1[(w_)?NumericQ, (r_)?NumericQ] := 
      NIntegrate[f1[w, r], {t, 0.001, r}];
Z2[(w_)?NumericQ, (r_)?NumericQ] := NIntegrate[f2[w, r], {t, boundary1, r}]

Obviously it should be

Z1[(w_)?NumericQ, (r_)?NumericQ] := 
          NIntegrate[f1[w, t], {t, 0.001, r}];
Z2[(w_)?NumericQ, (r_)?NumericQ] := NIntegrate[f2[w, t], {t, boundary1, r}]

After correction we have same picture for NDSolve and ParametricNDSolve (dashed lines) Figure 4

Improve this answer
$\endgroup$
8
  • $\begingroup$ Thank you for this. Can you see some obvious mistake that I have made when defining the functions using the solutions from ParametricNDSolve? I don't know another way of defining them for them to work as when I do so from NDSolve. $\endgroup$
    – juv95
    Commented Jan 29, 2023 at 19:41
  • 1
    $\begingroup$ More precisely, x2 in two definitions is not the same since Z1[1.48323, boundary1]=-0.431163, while Z[boundary1]=-0.230469. $\endgroup$
    – Alex Trounev
    Commented Jan 30, 2023 at 2:06
  • $\begingroup$ Thanks for the additional help. Since f1 are not different, all the issues have to be due to the integrating of it. Somehow in the second case(the ParametricNDSolve one), using NIntegrate on f1 is doing something different. I will try some stuff and post an update soon. $\endgroup$
    – juv95
    Commented Jan 30, 2023 at 5:34
  • 1
    $\begingroup$ @juv95 There is a typo in definition Z1[w, r] since you used under NIntegrate function f1[w,r] instead of f1[w,t]. $\endgroup$
    – Alex Trounev
    Commented Jan 30, 2023 at 6:19
  • $\begingroup$ Ahh, that was so careless of me... thank you for pointing it out. Everything is working out well now, except I had to make a change on the definition of Z1. Initially, I had it defined as NIntegrate[f1[w,t],{t,0.001,r}], and that was not working. When I switched it to NIntegrate[f1[w,r],{r,0.001,t}] instead, it works great. I am not sure why that is the case, but I'll mention it here in case it helps anyone in the future. $\endgroup$
    – juv95
    Commented Jan 30, 2023 at 7:34

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.