In this code, I am trying to solve a system of first-order ODE with corresponding boundary conditions. But it is showing some error which I am not able to rectify. Can anyone help me in rectifying the errors of this code. I will be very thankful to you. If any doubt is there please feel free to ask.
Thanking You in Advance
Constants
\[Beta] = 0.001;
\[Epsilon] = 0.1;
\[Alpha] = 1;
\[Zeta] = 0;
H = 1;
M = 1;
k = 1.5;
z = 0.5*\[Pi];
Subscript[\[Phi], 1] = 0.01;
Subscript[\[Phi], 2] = 0.01;
\[Rho]w = 997.1;
\[Mu]w = 0.001004;
\[Sigma]w = 5.5*10^-6;
\[Rho]k = 783;
\[Mu]k = 0.00151;
\[Sigma]k = 5*10^-11;
\[Rho]fe = 5200;
\[Sigma]fe = 25000;
\[Rho]Mi = 5060;
\[Sigma]Mi = 2.09*10^4;
t = \[Pi]/3;
Subscript[A, 1] = 1/(1 - Subscript[\[Phi], 1])^2.5;
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])^2.5*(1 - Subscript[\[Phi], 1])^2.5);
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;
Using NDSolve commands I am entering all the system of first order differential equations with corresponding boundary conditions
n1 = NDSolve[{(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) + ((][1]][1]
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, u100[1] == 0, u110[1] == 0,
u200[-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]), u100[0] == u200[0],
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)}, {u100, u110,
u200, u210, f1, f2, f3, f4, f5, f6, f7, f8, g1, g2, g3, g4, g5,
g6, g7, g8, g9, g10, g11, g12}, {y, -1, 1},
Method -> {"Shooting",
"ImplicitSolver" -> {"Newton",
"StepControl" -> "LineSearch"}} ];
Expressions
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;
Plot Command for the expression
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}}], {y, -1, 1}, PlotRange -> All, AspectRatio -> 1,
AxesLabel -> {y, U}, PlotStyle -> {Green, Thick}, ImageSize -> 260,
LabelStyle -> {FontSize -> 14, FontFamily -> "Times", Black, Bold}]