0
$\begingroup$

I just want to ask in the graph of the code mentioned in the link, while plotting u from {-0.0001,y,0.0001}, the jump is coming because of the numerical error u1-u2. Can anyone please tell me, how to minimise that error? So that the graph will become continuous.

Thanking You in advance Your help will be highly appreciable

Solving a system of ODE using shooting technique with NDSolve command

Code

\[Beta] = 10^-3;
\[Epsilon] = 1/10;
\[Alpha] = 1;
\[Zeta] = 0;
H = 1;
M = 1;
k = 3/2;
z = \[Pi]/2;
Subscript[\[Phi], 1] = 10^-2;
Subscript[\[Phi], 2] = 10^-2;
\[Rho]w = 997 + 1/10;
\[Mu]w = 1004*10^-6;
\[Sigma]w = 11/2*10^-6;
\[Rho]k = 783;
\[Mu]k = 151 10^-5;
\[Sigma]k = 5*10^-11;
\[Rho]fe = 5200;
\[Sigma]fe = 25000;
\[Rho]Mi = 5060;
\[Sigma]Mi = 209*10^2;
t = \[Pi]/3;
Subscript[A, 1] = 1/(1 - Subscript[\[Phi], 1])^(5/2);
Subscript[A, 
   2] = (1 - 
     Subscript[\[Phi], 1]) + ((Subscript[\[Phi], 1]*\[Rho]fe)/\[Rho]w);
Subscript[A, 
   3] = ((\[Sigma]fe*(1 + (2*
           Subscript[\[Phi], 1]))) + (2*\[Sigma]fe*(1 - 
         Subscript[\[Phi], 1])))/((\[Sigma]fe*(1 - 
         Subscript[\[Phi], 1])) + (\[Sigma]w*(2 + 
         Subscript[\[Phi], 1])));
Subscript[B, 1] = 
  1/((1 - Subscript[\[Phi], 2])^(5/2)*(1 - Subscript[\[Phi], 1])^(5/
        2));
Subscript[B, 
   2] = (((1 - 
         Subscript[\[Phi], 
          1]) + ((Subscript[\[Phi], 1]*\[Rho]fe)/\[Rho]k))*(1 - 
       Subscript[\[Phi], 2])) + (Subscript[\[Phi], 
       2]*\[Rho]Mi)/\[Rho]k;
Subscript[B, 
   3] = ((\[Sigma]Mi*(1 + (2*Subscript[\[Phi], 2]))) + (2*\[Sigma]w*
       Subscript[A, 
        3]*(1 - Subscript[\[Phi], 2])))/((\[Sigma]Mi*(1 - 
         Subscript[\[Phi], 2])) + (\[Sigma]w*
       Subscript[A, 3]*(2 + Subscript[\[Phi], 2])));
\[Sigma] = \[Sigma]w/\[Sigma]k;
\[Mu] = \[Mu]w/\[Mu]k;
\[Rho] = \[Rho]w/\[Rho]k;
M1 = Sqrt[\[Mu]/\[Sigma]]*M;
H1 = Sqrt[\[Mu]/\[Sigma]]*H;
eq = {(Subscript[A, 1]/
        Subscript[A, 
         2]*(u100''[
          y] + (1/(y + k)*u100'[y]) - (1/(y + k)^2*
           u100[y]))) - ((Subscript[A, 3]*M^2)/Subscript[A, 2]*
       u100[y]) - (k/(y + k)*H^2/Subscript[A, 2]) == 
    0, (Subscript[A, 1]/
        Subscript[A, 
         2]*(u110''[
          y] + (1/(y + k)*u110'[y]) - (1/(y + k)^2*
           u110[y]))) - (((Subscript[A, 3]*M^2)/Subscript[A, 2] + 
         I*H^2)*u110[y]) - (k/(y + k)*H^2/Subscript[A, 2]) == 
    0, (Subscript[B, 1]/
        Subscript[B, 
         2]*(u200''[
          y] + (1/(y + k)*u200'[y]) - (1/(y + k)^2*
           u200[y]))) - ((Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*
          M1^2)/Subscript[B, 2]*
       u200[y]) - (k/(y + k)*(\[Rho]*H1^2)/Subscript[B, 2]) == 
    0, (Subscript[B, 1]/
        Subscript[B, 
         2]*(u210''[
          y] + (1/(y + k)*u210'[y]) - (1/(y + k)^2*
           u210[y]))) - (((Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*
            M1^2)/Subscript[B, 2] + I*H1^2)*
       u210[y]) - (k/(y + k)*(\[Rho]*H1^2)/Subscript[B, 2]) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(f1''[y] + (1/(y + k)*f1'[y]))) - (f1[
        y]*((Subscript[A, 1]/
            Subscript[A, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
           Subscript[A, 2]))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(f2''[y] + (1/(y + k)*f2'[y]))) - (f2[
        y]*((Subscript[A, 1]/
            Subscript[A, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
           Subscript[A, 2]))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(f3''[y] + (1/(y + k)*f3'[y]))) - (f3[
        y]*((Subscript[A, 1]/
            Subscript[A, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
            Subscript[A, 2] + (I*H^2)))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(f4''[y] + (1/(y + k)*f4'[y]))) - (f4[
        y]*((Subscript[A, 1]/
            Subscript[A, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
            Subscript[A, 2] + (I*H^2)))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(f5''[y] + (1/(y + k)*f5'[y]))) - (f5[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]*
             Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(f6''[y] + (1/(y + k)*f6'[y]))) - (f6[
        y]*((Subscript[B, 1]/
            
            Subscript[B, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]*
             Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(f7''[y] + (1/(y + k)*f7'[y]))) - (f7[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]*
              Subscript[A, 3]*\[Sigma]*M1^2)/
            Subscript[B, 2] + (I*H1^2)))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(f8''[y] + (1/(y + k)*f8'[y]))) - (f8[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]*
              Subscript[A, 3]*\[Sigma]*M1^2)/
            Subscript[B, 2] + (I*H1^2)))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(g1''[y] + (1/(y + k)*g1'[y]))) - (g1[
        y]*((Subscript[A, 1]/
            Subscript[A, 2]*1/(y + k)^2) + (Subscript[A, 3]*M^2)/
          Subscript[A, 2])) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(g2''[y] + (1/(y + k)*g2'[y]))) - (g2[
        y]*((Subscript[A, 1]/
            Subscript[A, 
             2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[A, 3]*
             M^2)/Subscript[A, 2]))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(g3''[y] + (1/(y + k)*g3'[y]))) - (g3[
        y]*((Subscript[A, 1]/
            Subscript[A, 2]*(1/(y + k)^2 + 
             4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
           Subscript[A, 2]))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(g4''[y] + (1/(y + k)*g4'[y]))) - (g4[
        y]*((Subscript[A, 1]/
            Subscript[A, 2]*1/(y + k)^2) + ((Subscript[A, 3]*M^2)/
            Subscript[A, 2] + (I*H^2)))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(g5''[y] + (1/(y + k)*g5'[y]))) - (g5[
        y]*((Subscript[A, 1]/
            Subscript[A, 2]*(1/(y + k)^2 + 
             4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
            Subscript[A, 2] + (I*H^2)))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(g6''[y] + (1/(y + k)*g6'[y]))) - (g6[
        y]*((Subscript[A, 1]/
            Subscript[A, 2]*(1/(y + k)^2 + 
             4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
            Subscript[A, 2] + (I*H^2)))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(g7''[y] + (1/(y + k)*g7'[y]))) - (g7[
        y]*((Subscript[B, 1]/
            Subscript[B, 2]*(1/(y + k)^2)) + ((Subscript[B, 3]*
             Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(g8''[y] + (1/(y + k)*g8'[y]))) - (g8[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[B, 3]*
             Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(g9''[y] + (1/(y + k)*g9'[y]))) - (g9[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[B, 3]*
             Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(g10''[y] + (1/(y + k)*g10'[y]))) - (g10[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*1/(y + k)^2) + (((Subscript[B, 3]*
               Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I*
             H1^2)))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(g11''[y] + (1/(y + k)*g11'[y]))) - (g11[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + (((Subscript[B, 3]*
               Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I*
             H1^2)))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(g12''[y] + (1/(y + k)*g12'[y]))) - (g12[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + (((Subscript[B, 3]*
               Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I*
             H1^2)))) == 0};

U[-1] = ((u200[
        y] + \[Epsilon]*(f5[y]*Sin[\[Alpha]*z] + 
          f6[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g7[y] + 
          g8[y]*Sin[2*\[Alpha]*z] + 
          g9[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u210[
          y] + \[Epsilon]*(f7[y]*Sin[\[Alpha]*z] + 
            f8[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g10[y] + 
            g11[y]*Sin[2*\[Alpha]*z] + 
            g12[y]*Cos[2*\[Alpha]*z])))) /. y -> -1;

U[1] = ((u100[
        y] + \[Epsilon]*(f1[y]*Sin[\[Alpha]*z] + 
          f2[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g1[y] + 
          g2[y]*Sin[2*\[Alpha]*z] + 
          g3[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u110[
          y] + \[Epsilon]*(f3[y]*Sin[\[Alpha]*z] + 
            f4[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g4[y] + 
            g5[y]*Sin[2*\[Alpha]*z] + g6[y]*Cos[2*\[Alpha]*z])))) /. 
   y -> 1;

U1 = ((u200[
        y] + \[Epsilon]*(f5[y]*Sin[\[Alpha]*z] + 
          f6[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g7[y] + 
          g8[y]*Sin[2*\[Alpha]*z] + 
          g9[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u210[
          y] + \[Epsilon]*(f7[y]*Sin[\[Alpha]*z] + 
            f8[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g10[y] + 
            g11[y]*Sin[2*\[Alpha]*z] + 
            g12[y]*Cos[2*\[Alpha]*z])))) /. y -> 0;

U2 = ((u100[
        y] + \[Epsilon]*(f1[y]*Sin[\[Alpha]*z] + 
          f2[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g1[y] + 
          g2[y]*Sin[2*\[Alpha]*z] + 
          g3[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u110[
          y] + \[Epsilon]*(f3[y]*Sin[\[Alpha]*z] + 
            f4[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g4[y] + 
            g5[y]*Sin[2*\[Alpha]*z] + g6[y]*Cos[2*\[Alpha]*z])))) /. 
   y -> 0;

bc = {U[1] == 0, u110[1] == 0, U[-1] == 0, u210[-1] == 0, f2[1] == 0, 
   f4[1] == 0, f1[1] == -u100'[1], f3[1] == -u110'[1], 
   f5[-1] == -(Cos[\[Zeta]]*u200'[-1]), 
   f6[-1] == -(Sin[\[Zeta]]*u200'[-1]), 
   f7[-1] == -(Cos[\[Zeta]]*u210'[-1]), 
   f8[-1] == -(Sin[\[Zeta]]*u210'[-1]), 
   g1[1] == -(1/4*u100''[1]) - (1/2*f1'[1]), g2[1] == -(1/2*f2'[1]), 
   g3[1] == (1/4*u100''[1]) + (1/2*f1'[1]), 
   g4[1] == -(1/4*u110''[1]) - (1/2*f3'[1]), g5[1] == -(1/2*f4'[1]), 
   g6[1] == (1/4*u110''[1]) + (1/2*f3'[1]), 
   g7[-1] == -(Cos[\[Zeta]]/2*f5'[-1]) - (1/4*
       u200''[-1]) - (Sin[\[Zeta]]/2*f6'[-1]), 
   g8[-1] == -(Sin[\[Zeta]]/2*f5'[-1]) - (Sin[2*\[Zeta]]/4*
       u200''[-1]) - (Cos[\[Zeta]]/2*f6'[-1]), 
   g9[-1] == (Cos[\[Zeta]]/2*f5'[-1]) + (Cos[2*\[Zeta]]/4*
       u200''[-1]) - (Sin[\[Zeta]]/2*f6'[-1]), 
   g10[-1] == -(Cos[\[Zeta]]/2*f7'[-1]) - (1/4*
       u210''[-1]) - (Sin[\[Zeta]]/2*f8'[-1]), 
   g11[-1] == -(Sin[\[Zeta]]/2*f7'[-1]) - (Sin[2*\[Zeta]]/4*
       u210''[-1]) - (Cos[\[Zeta]]/2*f8'[-1]), 
   g12[-1] == (Cos[\[Zeta]]/2*f7'[-1]) + (Cos[2*\[Zeta]]/4*
       u210''[-1]) - (Sin[\[Zeta]]/2*f8'[-1]), U1 == U2, 
   u110[0] == 
    u210[0], ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(u100'[0] - 
       u100[0]/k) == (u200'[0] - 
      u200[0]/k), ((\[Mu]*Subscript[A, 1])/
       Subscript[B, 1])*(u110'[0] - u110[0]/k) == (u210'[0] - 
      u210[0]/k), f1[0] == f5[0], f2[0] == f6[0], f3[0] == f7[0], 
   f4[0] == 
    f8[0], ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f1'[0] - 
       f1[0]/k) == (f5'[0] - 
      f5[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f2'[0] - 
       f2[0]/k) == (f6'[0] - 
      f6[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f3'[0] - 
       f3[0]/k) == (f7'[0] - 
      f7[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f4'[0] - 
       f4[0]/k) == (f8'[0] - f8[0]/k), g1[0] == g7[0], g2[0] == g8[0],
    g3[0] == g9[0], g4[0] == g10[0], g5[0] == g11[0], 
   g6[0] == 
    g12[0], ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g1'[0] - 
       g1[0]/k) == (g7'[0] - 
      g7[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g2'[0] - 
       g2[0]/k) == (g8'[0] - 
      g8[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g3'[0] - 
       g3[0]/k) == (g9'[0] - 
      g9[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g4'[0] - 
       g4[0]/k) == (g10'[0] - 
      g10[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g5'[0] - 
       g5[0]/k) == (g11'[0] - 
      g11[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g6'[0] - 
       g6[0]/k) == (g12'[0] - g12[0]/k)};
var = {u100, u110, u200, u210, f1, f2, f3, f4, f5, f6, f7, f8, g1, g2,
   g3, g4, g5, g6, g7, g8, g9, g10, g11, g12}; coord = {y, -1, 
  1.2}; op = 
 Method -> {"Shooting", 
   "ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"}};

var0 = Table[var[[i]]@y, {i, Length[var]}]; var1 = D[var0, y]; var2 = 
 D[var1, y];

sol0 = Solve[eq, var2];

sol01 = sol0 /. y -> 1; sol02 = sol0 /. y -> -1;

bc1 = bc /. sol01[[1]] /. sol02[[1]];

Numerical solution

sol = NDSolve[{eq, bc1}, var, coord, op];

Boundary conditions test

bc1 /. sol[[1]]

Out[]= {False, False, False, True, False, False, False, False, \
True, True, True, False, False, False, False, False, False, False, \
True, True, True, True, True, True, False, False, False, False, \
False, False, False, False, False, False, False, False, False, False, \
False, False, False, False, False, False, False, False, False, False}

Velocity jump at y=0

(U1 - U2) /. sol[[1]]

Out[]= 1.62451*10^-7 + 1.48518*10^-10 I

The question is can we reduce the real part jump on one order to 10^-8?

Velocity jump visualization

With[{n1 = sol[[1]]}, 
  u101[y_] := (f1[y]*Sin[\[Alpha]*z] + f2[y]*Cos[\[Alpha]*z]) /. n1;
  u111[y_] := (f3[y]*Sin[\[Alpha]*z] + f4[y]*Cos[\[Alpha]*z]) /. n1;
  u201[y_] := (f5[y]*Sin[\[Alpha]*z] + f6[y]*Cos[\[Alpha]*z]) /. n1;
  u211[y_] := (f7[y]*Sin[\[Alpha]*z] + f8[y]*Cos[\[Alpha]*z]) /. n1;
  u102[y_] := (g1[y] + g2[y]*Sin[2*\[Alpha]*z] + 
      g3[y]*Cos[2*\[Alpha]*z]) /. n1;
  u112[y_] := (g4[y] + g5[y]*Sin[2*\[Alpha]*z] + 
      g6[y]*Cos[2*\[Alpha]*z]) /. n1;
  u202[y_] := (g7[y] + g8[y]*Sin[2*\[Alpha]*z] + 
      g9[y]*Cos[2*\[Alpha]*z]) /. n1;
  u212[y_] := (g10[y] + g11[y]*Sin[2*\[Alpha]*z] + 
      g12[y]*Cos[2*\[Alpha]*z]) /. n1;
  u10[y_] := (u100[y] + \[Epsilon]*u101[y] + \[Epsilon]^2*u102[y]) /. 
    n1;
  u11[y_] := (u110[y] + \[Epsilon]*u111[y] + \[Epsilon]^2*u112[y]) /. 
    n1;
  u20[y_] := (u200[y] + \[Epsilon]*u201[y] + \[Epsilon]^2*u202[y]) /. 
    n1;
  u21[y_] := (u210[y] + \[Epsilon]*u211[y] + \[Epsilon]^2*u212[y]) /. 
    n1];

p1 = Plot[
  Piecewise[{{Re[(u20[y] + (\[Beta]*(E^(I*t))*u21[y]))], -1 <= y <= 
      0}, {Re[(u10[y] + (\[Beta]*(E^(I*t))*u11[y]))], 
     0 <= y <= 1.2}}], {y, -10^-5, 10^-5}, PlotRange -> All, 
  AspectRatio -> 1, AxesLabel -> {"y", "U"}, ImageSize -> 260, 
  LabelStyle -> {FontSize -> 14, FontFamily -> "Times", Black, Bold}]

Figure 1

Visualization in 2D Figure 1

$\endgroup$
11
  • 1
    $\begingroup$ Please make your question self contained by including code and a picture $\endgroup$
    – MarcoB
    Commented Jun 9, 2023 at 13:19
  • $\begingroup$ Thank you, @MarcoB, for the suggestion, Actually the code is a little big, and already executed on Mathematica Stack Exchange. So, I think this is a better way to present the code. Kindly see by one simple click, please feel free to ask if you have any doubts, thank you so much. $\endgroup$ Commented Jun 9, 2023 at 16:59
  • $\begingroup$ @KomalGoyal In this model jump is about 1.62451*10^-7 only. Why do you need 10^-8? $\endgroup$ Commented Jun 11, 2023 at 2:31
  • $\begingroup$ @AlexTrounev, Because in the present case, if we plot 3-D function u(y,z) because of that jump that 3-D graph is discontinuous, but I want a continuous 3-D graph which will come when a numerical error is less than 10^-8 $\endgroup$ Commented Jun 11, 2023 at 3:45
  • $\begingroup$ @KomalGoyal The best what I got with NDSolve is u1-u2=3.83437*10^-8 + 3.25555*10^-11 I. $\endgroup$ Commented Jun 11, 2023 at 8:40

1 Answer 1

1
$\begingroup$

We can decrease error in 4.2 times by using NDSolve twice as follows

\[Beta] = 10^-3;
\[Epsilon] = 1/10;
\[Alpha] = 1;
\[Zeta] = 0;
H = 1;
M = 1;
k = 3/2;
Subscript[\[Phi], 1] = 10^-2;
Subscript[\[Phi], 2] = 10^-2;
\[Rho]w = 997 + 1/10;
\[Mu]w = 1004*10^-6;
\[Sigma]w = 11/2*10^-6;
\[Rho]k = 783;
\[Mu]k = 151 10^-5;
\[Sigma]k = 5*10^-11;
\[Rho]fe = 5200;
\[Sigma]fe = 25000;
\[Rho]Mi = 5060;
\[Sigma]Mi = 209*10^2;
t = \[Pi]/3;
Subscript[A, 1] = 1/(1 - Subscript[\[Phi], 1])^(5/2);
Subscript[A, 
   2] = (1 - 
     Subscript[\[Phi], 1]) + ((Subscript[\[Phi], 1]*\[Rho]fe)/\[Rho]w);
Subscript[A, 
   3] = ((\[Sigma]fe*(1 + (2*
           Subscript[\[Phi], 1]))) + (2*\[Sigma]fe*(1 - 
         Subscript[\[Phi], 1])))/((\[Sigma]fe*(1 - 
         Subscript[\[Phi], 1])) + (\[Sigma]w*(2 + 
         Subscript[\[Phi], 1])));
Subscript[B, 1] = 
  1/((1 - Subscript[\[Phi], 2])^(5/2)*(1 - Subscript[\[Phi], 1])^(5/
        2));
Subscript[B, 
   2] = (((1 - 
         Subscript[\[Phi], 
          1]) + ((Subscript[\[Phi], 1]*\[Rho]fe)/\[Rho]k))*(1 - 
       Subscript[\[Phi], 2])) + (Subscript[\[Phi], 
       2]*\[Rho]Mi)/\[Rho]k;
Subscript[B, 
   3] = ((\[Sigma]Mi*(1 + (2*Subscript[\[Phi], 2]))) + (2*\[Sigma]w*
       Subscript[A, 
        3]*(1 - Subscript[\[Phi], 2])))/((\[Sigma]Mi*(1 - 
         Subscript[\[Phi], 2])) + (\[Sigma]w*
       Subscript[A, 3]*(2 + Subscript[\[Phi], 2])));
\[Sigma] = \[Sigma]w/\[Sigma]k;
\[Mu] = \[Mu]w/\[Mu]k;
\[Rho] = \[Rho]w/\[Rho]k;
M1 = Sqrt[\[Mu]/\[Sigma]]*M;
H1 = Sqrt[\[Mu]/\[Sigma]]*H;

eq = {(Subscript[A, 1]/
        Subscript[A, 
         2]*(u100''[
          y] + (1/(y + k)*u100'[y]) - (1/(y + k)^2*
           u100[y]))) - ((Subscript[A, 3]*M^2)/Subscript[A, 2]*
       u100[y]) - (k/(y + k)*H^2/Subscript[A, 2]) == 
    0, (Subscript[A, 1]/
        Subscript[A, 
         2]*(u110''[
          y] + (1/(y + k)*u110'[y]) - (1/(y + k)^2*
           u110[y]))) - (((Subscript[A, 3]*M^2)/Subscript[A, 2] + 
         I*H^2)*u110[y]) - (k/(y + k)*H^2/Subscript[A, 2]) == 
    0, (Subscript[B, 1]/
        Subscript[B, 
         2]*(u200''[
          y] + (1/(y + k)*u200'[y]) - (1/(y + k)^2*
           u200[y]))) - ((Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*
          M1^2)/Subscript[B, 2]*
       u200[y]) - (k/(y + k)*(\[Rho]*H1^2)/Subscript[B, 2]) == 
    0, (Subscript[B, 1]/
        Subscript[B, 
         2]*(u210''[
          y] + (1/(y + k)*u210'[y]) - (1/(y + k)^2*
           u210[y]))) - (((Subscript[B, 3]*Subscript[A, 3]*\[Sigma]*
            M1^2)/Subscript[B, 2] + I*H1^2)*
       u210[y]) - (k/(y + k)*(\[Rho]*H1^2)/Subscript[B, 2]) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(f1''[y] + (1/(y + k)*f1'[y]))) - (f1[
        y]*((Subscript[A, 1]/
            Subscript[A, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
           Subscript[A, 2]))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(f2''[y] + (1/(y + k)*f2'[y]))) - (f2[
        y]*((Subscript[A, 1]/
            Subscript[A, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
           Subscript[A, 2]))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(f3''[y] + (1/(y + k)*f3'[y]))) - (f3[
        y]*((Subscript[A, 1]/
            Subscript[A, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
            Subscript[A, 2] + (I*H^2)))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(f4''[y] + (1/(y + k)*f4'[y]))) - (f4[
        y]*((Subscript[A, 1]/
            Subscript[A, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
            Subscript[A, 2] + (I*H^2)))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(f5''[y] + (1/(y + k)*f5'[y]))) - (f5[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]*
             Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(f6''[y] + (1/(y + k)*f6'[y]))) - (f6[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]*
             Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(f7''[y] + (1/(y + k)*f7'[y]))) - (f7[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]*
              Subscript[A, 3]*\[Sigma]*M1^2)/
            Subscript[B, 2] + (I*H1^2)))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(f8''[y] + (1/(y + k)*f8'[y]))) - (f8[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + \[Alpha]^2)) + ((Subscript[B, 3]*
              Subscript[A, 3]*\[Sigma]*M1^2)/
            Subscript[B, 2] + (I*H1^2)))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(g1''[y] + (1/(y + k)*g1'[y]))) - (g1[
        y]*((Subscript[A, 1]/
            Subscript[A, 2]*1/(y + k)^2) + (Subscript[A, 3]*M^2)/
          Subscript[A, 2])) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(g2''[y] + (1/(y + k)*g2'[y]))) - (g2[
        y]*((Subscript[A, 1]/
            Subscript[A, 
             2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[A, 3]*
             M^2)/Subscript[A, 2]))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(g3''[y] + (1/(y + k)*g3'[y]))) - (g3[
        y]*((Subscript[A, 1]/
            Subscript[A, 2]*(1/(y + k)^2 + 
             4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
           Subscript[A, 2]))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(g4''[y] + (1/(y + k)*g4'[y]))) - (g4[
        y]*((Subscript[A, 1]/
            Subscript[A, 2]*1/(y + k)^2) + ((Subscript[A, 3]*M^2)/
            Subscript[A, 2] + (I*H^2)))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(g5''[y] + (1/(y + k)*g5'[y]))) - (g5[
        y]*((Subscript[A, 1]/
            Subscript[A, 2]*(1/(y + k)^2 + 
             4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
            Subscript[A, 2] + (I*H^2)))) == 
    0, (Subscript[A, 1]/
        Subscript[A, 2]*(g6''[y] + (1/(y + k)*g6'[y]))) - (g6[
        y]*((Subscript[A, 1]/
            Subscript[A, 2]*(1/(y + k)^2 + 
             4*\[Alpha]^2)) + ((Subscript[A, 3]*M^2)/
            Subscript[A, 2] + (I*H^2)))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(g7''[y] + (1/(y + k)*g7'[y]))) - (g7[
        y]*((Subscript[B, 1]/
            Subscript[B, 2]*(1/(y + k)^2)) + ((Subscript[B, 3]*
             Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(g8''[y] + (1/(y + k)*g8'[y]))) - (g8[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[B, 3]*
             Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(g9''[y] + (1/(y + k)*g9'[y]))) - (g9[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + ((Subscript[B, 3]*
             Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(g10''[y] + (1/(y + k)*g10'[y]))) - (g10[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*1/(y + k)^2) + (((Subscript[B, 3]*
               Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I*
             H1^2)))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(g11''[y] + (1/(y + k)*g11'[y]))) - (g11[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + (((Subscript[B, 3]*
               Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I*
             H1^2)))) == 
    0, (Subscript[B, 1]/
        Subscript[B, 2]*(g12''[y] + (1/(y + k)*g12'[y]))) - (g12[
        y]*((Subscript[B, 1]/
            Subscript[B, 
             2]*(1/(y + k)^2 + (4*\[Alpha]^2))) + (((Subscript[B, 3]*
               Subscript[A, 3]*\[Sigma]*M1^2)/Subscript[B, 2]) + (I*
             H1^2)))) == 0};

U[-1] = ((u200[
        y] + \[Epsilon]*(f5[y]*Sin[\[Alpha]*z] + 
          f6[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g7[y] + 
          g8[y]*Sin[2*\[Alpha]*z] + 
          g9[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u210[
          y] + \[Epsilon]*(f7[y]*Sin[\[Alpha]*z] + 
            f8[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g10[y] + 
            g11[y]*Sin[2*\[Alpha]*z] + 
            g12[y]*Cos[2*\[Alpha]*z])))) /. y -> -1;

U[1] = ((u100[
        y] + \[Epsilon]*(f1[y]*Sin[\[Alpha]*z] + 
          f2[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g1[y] + 
          g2[y]*Sin[2*\[Alpha]*z] + 
          g3[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u110[
          y] + \[Epsilon]*(f3[y]*Sin[\[Alpha]*z] + 
            f4[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g4[y] + 
            g5[y]*Sin[2*\[Alpha]*z] + g6[y]*Cos[2*\[Alpha]*z])))) /. 
   y -> 1;

U1 = ((u200[
        y] + \[Epsilon]*(f5[y]*Sin[\[Alpha]*z] + 
          f6[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g7[y] + 
          g8[y]*Sin[2*\[Alpha]*z] + 
          g9[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u210[
          y] + \[Epsilon]*(f7[y]*Sin[\[Alpha]*z] + 
            f8[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g10[y] + 
            g11[y]*Sin[2*\[Alpha]*z] + 
            g12[y]*Cos[2*\[Alpha]*z])))) /. y -> 0;

U2 = ((u100[
        y] + \[Epsilon]*(f1[y]*Sin[\[Alpha]*z] + 
          f2[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g1[y] + 
          g2[y]*Sin[2*\[Alpha]*z] + 
          g3[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u110[
          y] + \[Epsilon]*(f3[y]*Sin[\[Alpha]*z] + 
            f4[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g4[y] + 
            g5[y]*Sin[2*\[Alpha]*z] + g6[y]*Cos[2*\[Alpha]*z])))) /. 
   y -> 0;

bc = {U[1] == 0, u110[1] == 0, U[-1] == 0, u210[-1] == 0, f2[1] == 0, 
   f4[1] == 0, f1[1] == -u100'[1], f3[1] == -u110'[1], 
   f5[-1] == -(Cos[\[Zeta]]*u200'[-1]), 
   f6[-1] == -(Sin[\[Zeta]]*u200'[-1]), 
   f7[-1] == -(Cos[\[Zeta]]*u210'[-1]), 
   f8[-1] == -(Sin[\[Zeta]]*u210'[-1]), 
   g1[1] == -(1/4*u100''[1]) - (1/2*f1'[1]), g2[1] == -(1/2*f2'[1]), 
   g3[1] == (1/4*u100''[1]) + (1/2*f1'[1]), 
   g4[1] == -(1/4*u110''[1]) - (1/2*f3'[1]), g5[1] == -(1/2*f4'[1]), 
   g6[1] == (1/4*u110''[1]) + (1/2*f3'[1]), 
   g7[-1] == -(Cos[\[Zeta]]/2*f5'[-1]) - (1/4*
       u200''[-1]) - (Sin[\[Zeta]]/2*f6'[-1]), 
   g8[-1] == -(Sin[\[Zeta]]/2*f5'[-1]) - (Sin[2*\[Zeta]]/4*
       u200''[-1]) - (Cos[\[Zeta]]/2*f6'[-1]), 
   g9[-1] == (Cos[\[Zeta]]/2*f5'[-1]) + (Cos[2*\[Zeta]]/4*
       u200''[-1]) - (Sin[\[Zeta]]/2*f6'[-1]), 
   g10[-1] == -(Cos[\[Zeta]]/2*f7'[-1]) - (1/4*
       u210''[-1]) - (Sin[\[Zeta]]/2*f8'[-1]), 
   g11[-1] == -(Sin[\[Zeta]]/2*f7'[-1]) - (Sin[2*\[Zeta]]/4*
       u210''[-1]) - (Cos[\[Zeta]]/2*f8'[-1]), 
   g12[-1] == (Cos[\[Zeta]]/2*f7'[-1]) + (Cos[2*\[Zeta]]/4*
       u210''[-1]) - (Sin[\[Zeta]]/2*f8'[-1]), U1 == U2, 
   u110[0] == 
    u210[0], ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(u100'[0] - 
       u100[0]/k) == (u200'[0] - 
      u200[0]/k), ((\[Mu]*Subscript[A, 1])/
       Subscript[B, 1])*(u110'[0] - u110[0]/k) == (u210'[0] - 
      u210[0]/k), f1[0] == f5[0], f2[0] == f6[0], f3[0] == f7[0], 
   f4[0] == 
    f8[0], ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f1'[0] - 
       f1[0]/k) == (f5'[0] - 
      f5[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f2'[0] - 
       f2[0]/k) == (f6'[0] - 
      f6[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f3'[0] - 
       f3[0]/k) == (f7'[0] - 
      f7[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(f4'[0] - 
       f4[0]/k) == (f8'[0] - f8[0]/k), g1[0] == g7[0], g2[0] == g8[0],
    g3[0] == g9[0], g4[0] == g10[0], g5[0] == g11[0], 
   g6[0] == 
    g12[0], ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g1'[0] - 
       g1[0]/k) == (g7'[0] - 
      g7[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g2'[0] - 
       g2[0]/k) == (g8'[0] - 
      g8[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g3'[0] - 
       g3[0]/k) == (g9'[0] - 
      g9[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g4'[0] - 
       g4[0]/k) == (g10'[0] - 
      g10[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g5'[0] - 
       g5[0]/k) == (g11'[0] - 
      g11[0]/k), ((\[Mu]*Subscript[A, 1])/Subscript[B, 1])*(g6'[0] - 
       g6[0]/k) == (g12'[0] - g12[0]/k)};

var = {u100, u110, u200, u210, f1, f2, f3, f4, f5, f6, f7, f8, g1, g2,
   g3, g4, g5, g6, g7, g8, g9, g10, g11, g12}; coord = {y, -1, 
  1.2}; op = 
 Method -> {"Shooting", 
   "ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"}};



var0 = Table[var[[i]]@y, {i, Length[var]}]; var1 = D[var0, y]; var2 = 
 D[var1, y];

sol0 = Solve[eq, var2];

sol01 = sol0 /. y -> 1; sol02 = sol0 /. y -> -1;

bc1 = bc /. sol01[[1]] /. sol02[[1]]; 

First numerical solution

sol = NDSolve[{eq, bc1}, var, coord, op];

With this solution we have

(U1 - U2) /. sol[[1]]

(*Out[]= 1.62451*10^-7 + 1.48518*10^-10 I*)

Second numerical solution

ic0 = var0 /. sol[[1]] /. y -> -1; ic01 = 
 var1 /. sol[[1]] /. y -> -1; ic1 = 
 Join[Table[var0[[i]] == ic0[[i]] /. y -> -1, {i, Length[var0]}], 
  Table[var1[[i]] == ic01[[i]] /. y -> -1, {i, Length[var1]}]];

 sol1 = 
 NDSolve[{eq, ic1}, var, coord, 
  Method -> {"DiscontinuityProcessing" -> False, 
    "TimeIntegration" -> "Extrapolation"}, 
  StartingStepSize -> 2 10^-13]; 

With this solution we compute

(U1 - U2) /. sol1[[1]]

Out[]= 3.83437*10^-8 + 3.25555*10^-11 I

Visualization

With[{n1 = sol[[1]]}, 
  u101[y_] := (f1[y]*Sin[\[Alpha]*z] + f2[y]*Cos[\[Alpha]*z]) /. n1;
  u111[y_] := (f3[y]*Sin[\[Alpha]*z] + f4[y]*Cos[\[Alpha]*z]) /. n1;
  u201[y_] := (f5[y]*Sin[\[Alpha]*z] + f6[y]*Cos[\[Alpha]*z]) /. n1;
  u211[y_] := (f7[y]*Sin[\[Alpha]*z] + f8[y]*Cos[\[Alpha]*z]) /. n1;
  u102[y_] := (g1[y] + g2[y]*Sin[2*\[Alpha]*z] + 
      g3[y]*Cos[2*\[Alpha]*z]) /. n1;
  u112[y_] := (g4[y] + g5[y]*Sin[2*\[Alpha]*z] + 
      g6[y]*Cos[2*\[Alpha]*z]) /. n1;
  u202[y_] := (g7[y] + g8[y]*Sin[2*\[Alpha]*z] + 
      g9[y]*Cos[2*\[Alpha]*z]) /. n1;
  u212[y_] := (g10[y] + g11[y]*Sin[2*\[Alpha]*z] + 
      g12[y]*Cos[2*\[Alpha]*z]) /. n1;
  u10[y_] := (u100[y] + \[Epsilon]*u101[y] + \[Epsilon]^2*u102[y]) /. 
    n1;
  u11[y_] := (u110[y] + \[Epsilon]*u111[y] + \[Epsilon]^2*u112[y]) /. 
    n1;
  u20[y_] := (u200[y] + \[Epsilon]*u201[y] + \[Epsilon]^2*u202[y]) /. 
    n1;
  u21[y_] := (u210[y] + \[Epsilon]*u211[y] + \[Epsilon]^2*u212[y]) /. 
    n1];

p1 = Plot[
  Piecewise[{{Re[(u20[y] + (\[Beta]*(E^(I*t))*u21[y]))], -1 <= y <= 
      0}, {Re[(u10[y] + (\[Beta]*(E^(I*t))*u11[y]))], 
     0 <= y <= 1.2}}], {y, -10^-5, 10^-5}, PlotRange -> All, 
  AspectRatio -> 1, AxesLabel -> {"y", "U"}, ImageSize -> 260, 
  LabelStyle -> {FontSize -> 14, FontFamily -> "Times", Black, Bold}, 
  Frame -> True]

With[{n1 = sol1[[1]]}, 
  u101[y_] := (f1[y]*Sin[\[Alpha]*z] + f2[y]*Cos[\[Alpha]*z]) /. n1;
  u111[y_] := (f3[y]*Sin[\[Alpha]*z] + f4[y]*Cos[\[Alpha]*z]) /. n1;
  u201[y_] := (f5[y]*Sin[\[Alpha]*z] + f6[y]*Cos[\[Alpha]*z]) /. n1;
  u211[y_] := (f7[y]*Sin[\[Alpha]*z] + f8[y]*Cos[\[Alpha]*z]) /. n1;
  u102[y_] := (g1[y] + g2[y]*Sin[2*\[Alpha]*z] + 
      g3[y]*Cos[2*\[Alpha]*z]) /. n1;
  u112[y_] := (g4[y] + g5[y]*Sin[2*\[Alpha]*z] + 
      g6[y]*Cos[2*\[Alpha]*z]) /. n1;
  u202[y_] := (g7[y] + g8[y]*Sin[2*\[Alpha]*z] + 
      g9[y]*Cos[2*\[Alpha]*z]) /. n1;
  u212[y_] := (g10[y] + g11[y]*Sin[2*\[Alpha]*z] + 
      g12[y]*Cos[2*\[Alpha]*z]) /. n1;
  u10[y_] := (u100[y] + \[Epsilon]*u101[y] + \[Epsilon]^2*u102[y]) /. 
    n1;
  u11[y_] := (u110[y] + \[Epsilon]*u111[y] + \[Epsilon]^2*u112[y]) /. 
    n1;
  u20[y_] := (u200[y] + \[Epsilon]*u201[y] + \[Epsilon]^2*u202[y]) /. 
    n1;
  u21[y_] := (u210[y] + \[Epsilon]*u211[y] + \[Epsilon]^2*u212[y]) /. 
    n1];

p2 = Plot[
  Piecewise[{{Re[(u20[y] + (\[Beta]*(E^(I*t))*u21[y]))], -1 <= y <= 
      0}, {Re[(u10[y] + (\[Beta]*(E^(I*t))*u11[y]))], 
     0 <= y <= 1.2}}], {y, -10^-5, 10^-5}, PlotRange -> All, 
  AspectRatio -> 1, AxesLabel -> {"y", "U"}, ImageSize -> 260, 
  LabelStyle -> {FontSize -> 14, FontFamily -> "Times", Black, Bold}, 
  PlotStyle -> {Red, Dashed}]

Show\[p1, p2\];][1]][1]

Figure 1

To plot velocity as a function y, z we use

f[y_, z_] := 
  If[y < 0, ((u200[
         y] + \[Epsilon]*(f5[y]*Sin[\[Alpha]*z] + 
           f6[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g7[y] + 
           g8[y]*Sin[2*\[Alpha]*z] + 
           g9[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u210[
           y] + \[Epsilon]*(f7[y]*Sin[\[Alpha]*z] + 
             f8[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g10[y] + 
             g11[y]*Sin[2*\[Alpha]*z] + g12[y]*Cos[2*\[Alpha]*z])))),
    ((u100[
         y] + \[Epsilon]*(f1[y]*Sin[\[Alpha]*z] + 
           f2[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g1[y] + 
           g2[y]*Sin[2*\[Alpha]*z] + 
           g3[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u110[
           y] + \[Epsilon]*(f3[y]*Sin[\[Alpha]*z] + 
             f4[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g4[y] + 
             g5[y]*Sin[2*\[Alpha]*z] + 
             g6[y]*Cos[2*\[Alpha]*z]))))] /. 
   NDSolve[{eq, bc1} /. z -> z, var, coord, op][[1]];

p2 = Plot[
  Evaluate[Table[Re[f[y, z]], {z, 0, 2 Pi, Pi/16}]], {y, -1, 1}, 
  PlotRange -> All, AspectRatio -> 1, AxesLabel -> {"y", "U"}, 
  ImageSize -> 260, 
  LabelStyle -> {FontSize -> 14, FontFamily -> "Times", Black, Bold}, 
  PlotLegends -> Table[z, {z, 0, 2 Pi, Pi/16}]]

Figure 2

To extract data from p2 we use

Z = Table[z, {z, 0, 2 Pi, Pi/16}] // N;
lst = Table[{Z[[i]], 
    Extract[p2[[1, 1]][[1]][[1]][[2 + i]][[1]][[2]], 1]}, {i, 
    Length[Z]}];
YY = Table[lst[[i, 2]][[All, 1]], {i, Length[Z]}];
uu = Table[lst[[i, 2]][[All, 2]], {i, Length[Z]}];
dimY = Table[Length[YY[[i]]], {i, Length[Z]}];

u = Interpolation[
   Flatten[Table[{{YY[[i, j]], Z[[i]]}, uu[[i, j]]}, {i, 
      Length[Z]}, {j, Min[dimY]}], 1], InterpolationOrder -> 1];

Plot3D[u[y, z], {y, -1, 1}, {z, 0, 2 Pi}, Mesh -> None, 
 ColorFunction -> "Aquamarine", PlotTheme -> "Marketing", 
 AxesLabel -> {"y", "z", "U"}]

Figure 3

Second method

Y1 = YY[[1]]; uz = ConstantArray[0, {Length[Y1], Length[Z]}];

Do[uz[[All, i]] = 
   Table[If[
       y < 0, ((u200[
            y] + \[Epsilon]*(f5[y]*Sin[\[Alpha]*z] + 
              f6[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g7[y] + 
              g8[y]*Sin[2*\[Alpha]*z] + 
              g9[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u210[
              y] + \[Epsilon]*(f7[y]*Sin[\[Alpha]*z] + 
                f8[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g10[y] + 
                g11[y]*Sin[2*\[Alpha]*z] + 
                g12[y]*Cos[2*\[Alpha]*z])))),
       ((u100[
            y] + \[Epsilon]*(f1[y]*Sin[\[Alpha]*z] + 
              f2[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g1[y] + 
              g2[y]*Sin[2*\[Alpha]*z] + 
              g3[y]*Cos[2*\[Alpha]*z])) + (\[Beta]*(E^(I*t))*(u110[
              y] + \[Epsilon]*(f3[y]*Sin[\[Alpha]*z] + 
                f4[y]*Cos[\[Alpha]*z]) + \[Epsilon]^2*(g4[y] + 
                g5[y]*Sin[2*\[Alpha]*z] + 
                g6[y]*Cos[2*\[Alpha]*z]))))] /. z -> Z[[i]], {y, 
      Y1}] /. NDSolve[{eq, bc1} /. z -> Z[[i]], var, coord, 
      op][[1]];, {i, Length[Z]}]

u2 = Interpolation[
  Flatten[Table[{{Y1[[j]], Z[[i]]}, uz[[j, i]]}, {i, Length[Z]}, {j, 
     dimY1[[2]]}], 1]]

Plot3D[u2[y, z] // Re, {y, -1, 1}, {z, 0, Pi}, Mesh -> None, 
 ColorFunction -> "Aquamarine", PlotTheme -> "Marketing", 
 AxesLabel -> {"y", "z", "U"}]

Figure 3

$\endgroup$
8
  • $\begingroup$ Thank You sir for time and consideration @AlexTrounev $\endgroup$ Commented Jun 11, 2023 at 12:41
  • $\begingroup$ @AlexTrounvev, can you please see the 3-D profiles for the same code, as mentioned above. $\endgroup$ Commented Jun 12, 2023 at 5:40
  • $\begingroup$ @KomalGoyal Do you mean Plot3D using z as coordinate? $\endgroup$ Commented Jun 12, 2023 at 6:09
  • $\begingroup$ @AlexTrounvev yes sir, I tried two times NDSolve also $\endgroup$ Commented Jun 12, 2023 at 7:11
  • 1
    $\begingroup$ @KomalGoyal Could you stop editing my answer? $\endgroup$ Commented Jun 12, 2023 at 17:28

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.