I am attempting to find the coefficients $c_1,c_2,c_3,c_4,c_5,c_6$ of the two functions $u(r) = c_1 r + c_2 (1/r) + c_3I_1(r/l_0) + c_4K_1(r/l_0)$ and $\phi(r) = c_5ln(r) + c_6 + (f/a)(\frac{du}{dr} + u(r)/r)$, where $I_1$ and $K_1$ are the modified Bessel function of the first and second kind, respectively, $l_0,f,a$ are constants. The boundary conditions to determined these coefficients are such that

  • $u(r_i)=u_1$
  • $u(r_o)=u_2$
  • $\phi(r_i)=\phi_1$
  • $\phi(r_o)=\phi_2$
  • $(\lambda+2\mu)\frac{du}{dr}+\lambda \frac{u}{r}+f(-\frac{d\phi}{dr})$ at $r_i$ = 0
  • $(\lambda+2\mu)\frac{du}{dr}+\lambda \frac{u}{r}+f(-\frac{d\phi}{dr})$ at $r_0$ = 0

The following is my Mathematica code:

Y = 139 10^9;
nu = 0.3;
lambda = (Y*nu)/((1 + nu)*(1 - 2*nu));
mu = 0.5*Y/(1 + nu);
f1 = 1 10^-6; f2 = 1 10^-6; f = f1 + 2*f2; kappa = 1 10^-9;
ls = 0.0 10^-6;
l0 = Sqrt[ls^2 + f^2/((lambda + mu)*kappa)];
dervI1[r_] := 1/2 (BesselI[0, r] + BesselI[2, r]);
dervI0[r_] := 1/2 (BesselI[-1, r] + BesselI[1, r]);
dervI2[r_] := 1/2 (BesselI[1, r] + BesselI[3, r]);
d2I1[r_] := 1/2 (dervI0[r] + dervI2[r]);
dervK1[r_] := -1/2 (BesselK[0, r] + BesselK[2, r]);
dervK0[r_] := -1/2 (BesselK[-1, r] + BesselK[1, r]);
dervK2[r_] := -1/2 (BesselK[1, r] + BesselK[3, r]);
d2K1[r_] := -1/2 (dervK0[r] + dervK2[r]);
u[r_] := c1 r + c2  1/r + c3 BesselI[1, r/l0] + c4 BesselK[1, r/l0];
du[r_] := c1 - c2 1/r^2 + 1/l0 c3 dervI1[r/l0] + 1/l0 c4 dervK1[r/l0];
ddu[r_] := -2 c2 1/r^3 + (1/l0)^2 c3 d2I1[r/l0] + (1/l0)^2 c4 d2K1[r/l0];
\[CurlyPhi][r_] := c5 Log[r] + c6 + f/kappa (du[r] + u[r]/r);
EE[r_] := -(c5 /r) - f/kappa (ddu[r] - u[r]/r^2 + du[r]/r);
\[Sigma]r[r_] := (lambda + 2*mu) du[r] + lambda u[r]/r + (f1 + 2*f2) EE[r];
d\[Sigma][r_] := D[\[Sigma]r[r], r];
eq1 = u[10 10^-6];
eq2 = u[20 10^-6];
eq3 = \[CurlyPhi][10 10^-6];
eq4 = \[CurlyPhi][20 10^-6];
eq5 = l0^2*d\[Sigma][r] /. (r -> 10 10^-6);
eq6 = l0^2*d\[Sigma][r] /. (r -> 20 10^-6);
uin = 0.045 10^-6;
uout = 0.05 10^-6;
\[CurlyPhi]in = 0;
\[CurlyPhi]out = 1;
Sol = NSolve[{eq1 == uin, eq2 == uout, eq3 == \[CurlyPhi]in, eq4 == \[CurlyPhi]out, eq5 == 0, eq6 == 0}, {c1, c2, c3, c4, c5,c6}];
c1 = Sol[[1]][[1]][[2]];
c2 = Sol[[1]][[2]][[2]];
c3 = Sol[[1]][[3]][[2]];
c4 = Sol[[1]][[4]][[2]];
c5 = Sol[[1]][[5]][[2]];
c6 = Sol[[1]][[6]][[2]];
(*root=FindRoot[{eq1==uin,eq2==uout,eq3==\[CurlyPhi]in,eq4== 
\[CurlyPhi]out,eq5==0,eq6==0},{{c1,0.0018069515755701755},{c2,2.\
8063598198085684*^-13},{c3,-6.70326128908629*^-14},{c4,-2.\
671175261048626*^-7},{c5,4.016694593136571},{c6,33.86980711923503}},\
AccuracyGoal->Infinity,PrecisionGoal->8]
c1=root[[1]][[2]];
c2=root[[2]][[2]];
c3=root[[3]][[2]];
c4=root[[4]][[2]];
c5=root[[5]][[2]];
c6=root[[6]][[2]];*)
Plot[u[r], {r, 10 10^-6, 20 10^-6}]
Plot[\[CurlyPhi][r], {r, 10 10^-6, 20 10^-6}]

For some other values of $l_0$ I can get the correct coefficients, but for the value in the above code, I could not. Should I use different function? I appreciate any help or suggestion. Thank you for reading.

Update

I have solved the problem by using NSolve and SetPrecision

Sol = NSolve[{SetPrecision[eq1 == uin, 500], SetPrecision[eq2 == uout, 500], SetPrecision[eq3 == \[CurlyPhi]in, 500], SetPrecision[eq4 == \[CurlyPhi]out, 500], SetPrecision[eq5 == 0, 500], SetPrecision[eq6 == 0, 500]}, {c1, c2, c3, c4, c5, c6}, Reals];

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.