2
$\begingroup$

My code is given below. I am getting each and every plot from the start. But I am unable to get the final plot in the code shown below.

T[B_] := NIntegrate[(Exp[-(B/(f0[z])^2)*t] - 1)/t, {t, 0, 1}];
T1[r_, b_, B_, W_, z_] := NIntegrate[t^((r^2*(f0[z])^2)/(b^2*(f[z])^2) - 1)*(1 -Exp[-B*t]/(f0[z])^2), {t, 0, 1}];
T2[r_, b_, B_, W_, z_] := NIntegrate[(Log[t])*(1 - Exp[-B*t]/(f0[z])^2)*t^((r^2*(f0[z])^2)/(b^2*(f[z])^2) - 1), {t, 0, 1}];
T3[r_, b1_, B_, W_, z_] := NIntegrate[t^((r^2*(f0[z])^2)/((b1)^2*(fs[z])^2) - 1)*(1 - Exp[-B*t]/(f0[z])^2), {t, 0, 1}];
T4[r_, b1_, B_, W_, z_] :=NIntegrate[(Log[t])*(1 - Exp[-B*t]/(f0[z])^2)*t^((r^2*(f0[z])^2)/((b1)^2*(fs[z])^2) - 1), {t, 0, 1}];
S[wp_, L_] := 1/2*Sqrt[\[Pi]/8]*wp/L^3*Exp[-1/(2*L^2) - 3/2];
ki[wp_, L_, w_, k_, c_] := (S[wp, L]*w)/(k*(0.1)^2*c^2);
s[r_, b_, b1_, B_, W_, G_, wp_, L_, w_, k0_, k_, ks0_, c_, z_] :=NDSolve[{f0''[z] + 1/f0[z]*(f0'[z])^2 == 1/(f0[z])^3 - 
  W^2*1/f0[z]*((f0[z])^2/B*(Exp[-B/(f0[z])^2] - 1 - T[B])), f''[z] + 1/f[z]*(f'[z])^2 == 1/(4*(f0[z])^3) - 1/(4*G)* W^2*((f0[z])^2/(f[z])^2*T1[r, b, B, W, z] + (f0[z])^4/(f[z])^4*T2[r, b, B, W, z])(**Exp[2*ki[wp,L,w,k,c]*z*k0*r^2]*)(*(S[
wp,L]^2*w^2*f[z]*r^4)/(2*G^2*(0.1)^2*c^2*Exp[-2*ki[wp,L,w,k,c]*z*
k0*r^2])*), fs''[z] + 1/fs[z]*(fs'[z])^2 == (k0^2*r^4)/((ks0)^2*(b1)^4*(fs[z])^3) - (k0^2*r^4*(f0[z])^2)/(ks0^2*b1^4*(fs[z])^3)*W^2*(T3[r, b1, B, W, z] + (r^2*(f0[z])^2)/(b1^2*(fs[z])^2)*T4[r, b1, B, W, z]), f0[0] == 1, f0'[0] == 0, f[0] == 1,f'[0] == 0, fs[0] == 1, fs'[0] == 0}, {f0, f, fs}, {z, 0, 3.0}, Method -> {"IndexReduction" -> Automatic}];
a[z_] = s[24*10^(-6), 24*10^(-6), 48*10^(-6), 1.5, 2.5, 3,0.306*10^14, 0.2, 10^14, 0.59*10^7, 1.18*10^7, 0.287*10^7, 3*10^8, z] // Quiet
F0G = Evaluate[f0[z] /. a[z]]
FG = Evaluate[f[z] /. a[z]]
FSG = Evaluate[fs[z] /. a[z]]
F0 = Evaluate[f0[0.212] /. a[z]]
F = Evaluate[f[0.212] /. a[z]]
Fs = Evaluate[fs[0.212] /. a[z]]
I1[B_, r_, b1_, w0_, wp_, w_, ws_, L_, k_, k0_, ks0_, ks1_, c_, z_] :=1/4*(wp^4*ws^2)/(c^4*w0^2)*((7.6*10^-3))^2*1/((ks1)^2 - (ks0)^2 -wp^2/c^2*(1 - Exp[-(B/(F0)^2)*Exp[-((b1*Fs)/(r*F0))^2]]))^2*((Fs/(F0*F))^2)*1/(FSG)^2*Exp[-2*ki[wp, w, L, k, c]*0.212*k0*r^2 - 1];
CI1[z_] = I1[1.5, 24*10^(-6), 48*10^(-6), 1.77*10^14, 0.306*10^14, 10^14,3.54*10^14, 0.2, 1.18*10^7, 0.59*10^7, 0.287*10^7, -0.59*10^7,3*(10^8), z]
Plot[CI1[z], {z, 0, 3}, PlotLegends -> "Expressions",PlotRange -> Full]
I2[B_, r_, b_, b1_, w0_, wp_, w_, ws_, L_, k_, k0_, ks0_, ks1_, c_,z_] := 1/4*(wp^4*ws^2)/(c^4*w0^2)*((7.6*10^-3))^2*1/((ks1)^2 - (ks0)^2 - wp^2/c^2*(1 Exp[-(B/(F0)^2)*Exp[-((b1*Fs)/(r*F0))^2]]))^2*1/((FG)^2*(F0G)^2)*Exp[-((b1*FSG)/(b*FG))^2 - ((b1*FSG)/(r*F0G))^2]*Exp[-2*ki[wp, w, L, k, c]*0.212*k0*r^2];
CI2[z_] = I2[1.5, 24*10^(-6), 24*10^(-6), 48*10^(-6), 1.77*10^14, 0.306*10^14,10^14, 3.54*10^14, 0.2, 1.18*10^7, 0.59*10^7, 0.287*10^7, -0.59*10^7, 3*(10^8), z]
Plot[CI2[z], {z, 0, 3}, PlotLegends -> "Expressions", PlotRange -> Full]
I3[B_, r_, b_, b1_, w0_, wp_, w_, ws_, L_, k_, k0_, ks0_, ks1_, c_, z_] := -2*(wp^4*ws^2)/(c^4*w0^2)*((7.6*10^-3))^2*1/((ks1)^2 - (ks0)^2 - wp^2/c^2*(1 - Exp[-(B/(F0)^2)*Exp[-((b1*Fs)/(r*F0))^2]]))^2*(Fs/ F0*F)*Exp[-(1/2) - 1/2*((b1*FSG)/(b*FG))^2 -1/2*((b1*FSG)/(r*F0G))^2 -ki[wp, w, L, k, c]*(0.212*k0*r^2 +z*k0*r^2)]*1/(F0G*FG*FSG)*Cos[z*r^2*(ks0 + ks1) - 0.212*k0*r^2*(ks0 + ks1)];
CI3[z_] =I3[1.5, 24*10^(-6), 24*10^(-6), 48*10^(-6), 1.77*10^14, 0.306*10^14,10^14, 3.54*10^14, 0.2, 1.18*10^7, 0.59*10^7,0.287*10^7, -0.59*10^7, 3*(10^8), z]
Plot[CI3[z], {z, 0, 3}, PlotLegends -> "Expressions",PlotRange -> Full]

All is going well up to this point. But the code for reflectivity given below doesn't give me any plot.

R[B_, W_, r_, b_, b1_, w0_, wp_, w_, ws_, L_, k_, k0_, ks0_, ks1_, c_,z_] := I1[B, r, b, b1, w0, wp, w, ws, L, k, k0, ks0, ks1, c, z] - I2[B, r, b, b1, w0, wp, w, ws, L, k, k0, ks0, ks1, c, z] -I3[B, r, b, b1, w0, wp, w, ws, L, k, k0, ks0, ks1, c, z];
Ref[z_] = R[1.5, 2.5, 24*10^(-6), 24*10^(-6), 48*10^(-6), 1.77*10^14,0.306*10^14, 10^14, 3.54*10^14, 0.2, 1.18*10^7, 0.59*10^7,0.287*10^7,-0.59*10^7, 3*(10^8), z]
Plot[Ref[z], {z, 0, 3}, PlotLegends -> "Expressions", PlotRange -> Full]
$\endgroup$
2
  • 4
    $\begingroup$ Could you do some basic troubleshooting on your own first? Where does the code go wrong? Try evaluating each result as you go, remove Quiet and pay attention to any warning, then evaluate the expression you are trying to plot to see if it returns real numbers. Finally, there are many Plot expressions in your code. Which one fails? $\endgroup$
    – MarcoB
    Feb 22 at 13:55
  • $\begingroup$ I am unable to get it on my own. Plot for R1 fails. $\endgroup$ Feb 23 at 6:05
4
$\begingroup$

You are not getting a plot because you have serious numerics problems in your code. The values of the functions FE, FS and R1 are often far beyond the range that can be represented by machine arithmetic. You need to re-express your code so it can make use of the Wolfram Languages arbitrary precision arithmetic. I am not going to do that for you. What I will do is re-express you code so it is simpler and makes the numeric problems more obvious.

s[wp_, w0_, B_, n_, z_] := 
  NDSolveValue[
    {(1 - (wp^2/w0^2)*Exp[-B/(f0[z])^2*n^(2*n)*Exp[-2*n]])*f0[z]*f0''[z] == 
      4/(f0[z])^2 - 
        2*(wp^2/w0^2)*(B/(f0[z])^2)*n^(2*n)*Exp[-2*n]*
          Exp[-(B/(f0[z])^2)*n^(2*n)*Exp[-2*n]]*(3481)^2, 
      f0[0] == 1, 
    f0'[0] == 0}, 
   f0, {z, 0, .01}];
F = s[0.79*10^15, 1.77*10^15, 5, 1, z]

interpF

Plot[F[z], {z, 0, 0.01}]

f_plot

Every thing is OK at this point.

ss[wp_, w0_, B_, n_, r_, b_, c_, k_, z_] :=
  NDSolveValue[
   {fe''[z] ==
    ((fe[z]*(3481)^2)/F[z]^2)*
     ((1/(k^2*r^2*F[z]^2)*(3 + ((r*F[z])/(b*fe[z]))^4)) - 
        2*wp^2/(3*k^2*0.1^2*c^2)*B/F[z]^2*n^(2*n)*Exp[-2*n]*
          Exp[-B/F[z]^2*n^(2*n)*Exp[-2*n]]),
   fe[0] == 1, fe'[0] == 0}, fe, {z, 0, .01}]
FE = ss[0.79*10^15, 1.77*10^15, 5, 1, 10*10^-6, 10*10^-6, 3*10^8, 1.18*10^7, z];
Plot[FE[z], {z, 0, 0.01}, PlotRange -> All]

fe_plot

At this point the numerics are already unsafe.

sss[wp_, w0_, ws_, B_, n_, r_, b_, b1_, c_, ks0_, z_] := 
  NDSolveValue[
    {fs''[z] == 
      (fs[z]*(3481)^2)/F[z]^2*
        (1/(ks0^2*r^2*F[z]^2)*(3 + ((r*F[z])/(b1*fs[z]))^4)) - 
          ws^2/(ks0^2*c^2)*2*(wp^2/ws^2*B/F[z]^2*n^(2*n)*Exp[-2*n])*
            Exp[-B/F[z]^2*n^(2*n)*Exp[-2*n]], 
    fs[0] == 1, fs'[0] == 0}, 
   fs, {z, 0, .01}]
FS = 
  sss[0.79*10^15, 1.77*10^15, 1.24*10^15, 5, 1, 10*10^-6, 
      10*10^-6, 10*10^-6, 3*10^8, 0.287*10^7, z]

At this point machine arithmetic is failing big time. Mathemaica is complaining bittely about it. FS is unusable for furth0er computation.

Another problem with your code is that η needs to be defined as a function of z. I would rewrite it as

Clear[η]
η[b1_, r_, n_, z_] := (b1*FS[z])/(r*F[z]) - (2*n)^(1/2)

then

η[10*10^-6, 10*10^-6, 1, z]

evaluates to something usable. R1 will also have to be rewritten.

R[wp_, w0_, w_, ws_, B_, n_, r_, b_, b1_, L_, c_, k_, ks0_, ks1_, n0_, N0_, z_] := 
  1/4*(wp^2/c^2)*(n0/N0)^2*(ws/w0)^2*
    1/
      (ks1^2 - ks0^2 - 
         wp^2/c^2*Exp[-B/F[z]^2((b1*FS[z])^2/(2*(r*F[z])^2))^n*
           Exp[-((b1^2*FS[z]^2)/(2*r^2*F[z]^2))]])^2*
    ((r*F[z])/(b*FE[z]))^(4*n)*
      (η[b1, r, n, z] + 2*n)^(1/2))^(8*n)/(2^(4*n)*FE[z]^2*F[z]^2)*
        Exp[-(η[b1, r, n, z] + (2*n)^(1/2))^2*((r*F[z])/(b*FE[z]))^2 - 
          (η[b1, r, n, z] + (2*n)^(1/2))^2 - 2*ki[wp, w, L, k, c]*z]

However, R1 remains unusable because of the numerics problems of its constituents, which is why you can't plot it.

$\endgroup$
2
  • $\begingroup$ Do you mean that there might be something wrong with the arguments which I give to my parameters? $\endgroup$ Feb 24 at 16:34
  • $\begingroup$ @amanbhatia. No, my point is that your can't use machine arithmetic to numerically solve your differential equations. $\endgroup$
    – m_goldberg
    Feb 26 at 16:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.