How to solve this set of integrated PDEs

Question

Could anyone help me solve set of pdes with embedded integration? The pdes below is solvable , but can not be solved when I tried different parameters such as B=10 *Exp[-r] or T=1000*Exp[-r] orn=100*Exp[-r] orn=100*(1-Exp[-r]).

ClearAll["Global`*"] ;
SetOptions[Plot, BaseStyle -> {FontSize -> 16, FontWeight -> Bold}, PlotTheme -> "Detailed", ImageSize -> Medium];
λ0 = 1.064; k0 = 2 π/λ0;b0 = 10^4; qr0 = (k0*b0)^-2;

equ = With[{ν = ν@z, qr = qr@z, qi = qi@z, n = 100, B = 10,TeV = 100, Z = 5/2, λ = 1, rMax = 10^5},
      With[{νe = 1.72*Z*n*λ^-1*TeV^(-3/2)},
      With[{ϵr =  1 - (n (1 - B ))/((1 - B )^2 + νe^2), ϵi = (nνe)/((1 - B )^2 + νe^2)},
      With[{ 
       K1 = Integrate[((1 - ϵr)*Cos[qi r^2] - ϵi*Sin[qi r^2]) Exp[-qr r^2] r, {r, 0, rMax}],
       K2 =Integrate[((1 - ϵr)*Sin[qi r^2] + ϵi* Cos[qi r^2]) Exp[-qr r^2] r, {r, 0, rMax}], 
       K3 = Integrate[ϵi*Exp[-2 qr*r^2] r, {r, 0, rMax}]},
       {ν D[qi, z] - qi D[ν, z] + ν (K1 (qi^2 - qr^2) - 2 K2 qi qr) == 0,
        ν D[qr, z] - qr D[ν, z] + ν (K2 (qi^2 - qr^2) + 2 K1 qi qr) == 0,
        -ν D[qr, z] + 2 qr D[ν, z] + 4 K3 qr^2 ν == 0}]]]];
ic = {ν[0] == 1, qr[0] == qr0, qi[0] == 0}
{νsol, qrsol, qisol} = NDSolveValue[{equ, ic}, {ν, qr, qi}, {z, 0, 10^2}];

zMax = 10^-1;
p1 = Plot[{νsol[k0*z]}, {z, 0, zMax}, PlotRange -> {{0, zMax}, {0, 1}}, FrameLabel -> {"z (μm)", "ν"}];
p2 = Plot[{qrsol[k0*z]^(-1/2)/k0}, {z, 0, zMax}, PlotRange -> {{0, zMax}, {0, qr0^(-1/2)/k0*2}},  FrameLabel -> {"z (μm)", "b (mm)"}];
Grid[{{p1, p2}}]

In this way, I revised the code using NIntegrate as below, but also failed!!!

ClearAll["Global`*"] ;
SetOptions[Plot, BaseStyle -> {FontSize -> 16, FontWeight -> Bold}, PlotTheme -> "Detailed", ImageSize -> Medium];

λ0 = 1.064; k0 = 2 π/λ0; b0 = 10^4; qr0 = (k0*b0)^-2;
n = 100; B = 10; TeV = 1000 *E^-r; Z = 5/2; λ = 1; rMax = 10^5;
νe = 1.72*Z*n*λ^-1*TeV ^(-3/2);
ϵr = 1 - (n (1 - B ))/((1 - B )^2 + νe^2);
ϵi = (n νe)/((1 - B )^2 + νe^2);

coreA[ν_, qr_, qi_, r_] := ((1 - ϵr)*Cos[qi r^2] - ϵi* Sin[qi r^2])Exp[-qr r^2] r
coreB[ν_, qr_, qi_, r_] := ((1 - ϵr)*Sin[qi r^2] + ϵi*Cos[qi r^2]) Exp[-qr r^2] r
coreC[qr_, r_] := ϵi*Exp[-2 qr*r^2] r

K1[ν_, qr_, qi_?NumericQ] := NIntegrate[coreA[ν, qr, qi, r], {r, 0, rMax}]
K2[ν_, qr_, qi_?NumericQ] := NIntegrate[coreB[ν, qr, qi, r], {r, 0, rMax}]
K3[qr_?NumericQ] := NIntegrate[coreC[qr, r], {r, 0, rMax}]

zMin = 0;
equ = With[{ν = ν@z, qr = qr@z, qi = qi@z},
     With[{K1 = K1[ν, qr, qi], K2 = K2[ν, qr, qi], K3 = K3[qr]},
  { D[ν, z] == -2 K1 ν qi - K2 ν qr^-1 (qi^2 - qr^2) -  4 K3 ν qr,
    D[qr, z] == -4 K1 qr qi - 2 K2 (qi^2 - qr^2) - 4 K3 qr^2,
    D[qi, z] == -K1 (3 qi^2 - qr^2) + K2 (-qi^3 qr^-1 + 3 qi qr) - 4 K3 qi qr}
]];
ic = {ν[zMin] == 1, qr[zMin] == qr0, qi[zMin] == 0}
{νsol, qrsol, qisol} =  NDSolveValue[{equ, ic}, {ν, qr, qi}, {z, zMin, 10^3}];

 zMax = 10^1;

 p1 = Plot[{νsol[k0*z]}, {z, 0, zMax}, PlotRange -> {{0, zMax}, {0, 1}}, FrameLabel -> {"z (μm)", "ν"}];
 p2 = Plot[{qrsol[k0*z]^(-1/2)/k0}, {z, 0, zMax}, PlotRange -> {{0, zMax}, {0, qr0^(-1/2)/k0*2}}, 
 FrameLabel -> {"z (μm)", "b (mm)"}];
 Grid[{{p1, p2}}]

Could any one help me? THx!

KL

Answer 1

We can use combination of FDM and iteration method as follows. First we compute test example to compare

ClearAll["Global`*"];
\[Lambda]0 = 1.064; k0 = 2 \[Pi]/\[Lambda]0; b0 = 
 10^4; qr0 = (k0*b0)^-2; equ = 
 With[{\[Nu] = \[Nu]@z, qr = qr@z, qi = qi@z, n = 100, B = 10, 
   TeV = 100, Z = 5/2, \[Lambda] = 1, rMax = 10^5}, 
  With[{\[Nu]e = 1.72*Z*n*\[Lambda]^-1*TeV^(-3/2)}, 
   With[{\[Epsilon]r = 
      1 - (n (1 - 
            B))/((1 - B)^2 + \[Nu]e^2), \[Epsilon]i = (n \[Nu]e)/((1 -
             B)^2 + \[Nu]e^2)}, 
    With[{K1 = 
       Integrate[((1 - \[Epsilon]r)*Cos[qi r^2] - \[Epsilon]i*
            Sin[qi r^2]) Exp[-qr r^2] r, {r, 0, rMax}], 
      K2 = Integrate[((1 - \[Epsilon]r)*Sin[qi r^2] + \[Epsilon]i*
            Cos[qi r^2]) Exp[-qr r^2] r, {r, 0, rMax}], 
      K3 = Integrate[\[Epsilon]i*Exp[-2 qr*r^2] r, {r, 0, 
         rMax}]}, {\[Nu] D[qi, z] - 
        qi D[\[Nu], z] + \[Nu] (K1 (qi^2 - qr^2) - 2 K2 qi qr) == 
       0, \[Nu] D[qr, z] - 
        qr D[\[Nu], z] + \[Nu] (K2 (qi^2 - qr^2) + 2 K1 qi qr) == 
       0, -\[Nu] D[qr, z] + 2 qr D[\[Nu], z] + 4 K3 qr^2 \[Nu] == 
       0}]]]];
ic = {\[Nu][0] == 1, qr[0] == qr0, qi[0] == 0};
{\[Nu]sol0, qrsol0, qisol0} = 
  NDSolveValue[{equ, ic}, {\[Nu], qr, qi}, {z, 0, 10}];

Now we combine FDM and iteration method in one code. Define parameters (same as above) and functions to compute with FDM

 n = 100; B = 10; TeV = 100; Z = 5/2; \[Lambda] = 1; rMax = 
 10^5; zMin = 0;
\[Nu]e = 1.72*Z*n*\[Lambda]^-1*TeV^(-3/2);
\[Epsilon]r = 1 - (n (1 - B))/((1 - B)^2 + \[Nu]e^2);
\[Epsilon]i = (n \[Nu]e)/((1 - B)^2 + \[Nu]e^2);

coreA[\[Nu]_, qr_, qi_, 
  r_] := ((1 - \[Epsilon]r)*Cos[qi r^2] - \[Epsilon]i*
     Sin[qi r^2]) Exp[-qr r^2] r
coreB[\[Nu]_, qr_, qi_, 
  r_] := ((1 - \[Epsilon]r)*Sin[qi r^2] + \[Epsilon]i*
     Cos[qi r^2]) Exp[-qr r^2] r
coreC[qr_, r_] := \[Epsilon]i*Exp[-2 qr*r^2] r

K1[\[Nu]_?NumericQ, qr_?NumericQ, qi_?NumericQ] := 
 NIntegrate[coreA[\[Nu], qr, qi, r], {r, 0, rMax}]
K2[\[Nu]_?NumericQ, qr_?NumericQ, qi_?NumericQ] := 
 NIntegrate[coreB[\[Nu], qr, qi, r], {r, 0, rMax}]
K3[qr_?NumericQ] := NIntegrate[coreC[qr, r], {r, 0, rMax}];

Initial data for iteration and loop to make 10 iterations

\[Nu]sol[0][z_] := 1; qrsol[0][z_] := qr0; qisol[0][z_] := 0;

Do[k1 = Interpolation[
   Table[{z, 
     K1[\[Nu]sol[i - 1][z], qrsol[i - 1][z], qisol[i - 1][z]]}, {z, 0,
      10, .1}]]; 
 k2 = Interpolation[
   Table[{z, 
     K2[\[Nu]sol[i - 1][z], qrsol[i - 1][z], qisol[i - 1][z]]}, {z, 0,
      10, .1}]]; 
 k3 = Interpolation[
   Table[{z, K3[qrsol[i - 1][z]]}, {z, 0, 10, .1}]]; {\[Nu]sol[i], 
   qrsol[i], qisol[i]} = 
  NDSolveValue[{Derivative[1][\[Nu]][z] == -2 k1[z] qi[z] \[Nu][z] - 
      4 k3[z] qr[z] \[Nu][z] - (k2[z] (qi[z]^2 - qr[z]^2) \[Nu][z])/
      qr[z], Derivative[1][qr][z] == -4 k1[z] qi[z] qr[z] - 
      4 k3[z] qr[z]^2 - 2 k2[z] (qi[z]^2 - qr[z]^2), 
    Derivative[1][qi][z] == -4 k3[z] qi[z] qr[z] + 
      k2[z] (-(qi[z]^3/qr[z]) + 3 qi[z] qr[z]) - 
      k1[z] (3 qi[z]^2 - qr[z]^2), \[Nu][zMin] == 1, qr[zMin] == qr0, 
    qi[zMin] == 0}, {\[Nu], qr, qi}, {z, zMin, 10}];, {i, 10}]

Finally we can compare numerical and exact solution

Table[{Plot[{\[Nu]sol[i][k0*z], \[Nu]sol0[k0*z]}, {z, 0, zMax}, 
   PlotRange -> All, FrameLabel -> {"z (\[Mu]m)", "\[Nu]"}, 
   PlotLabel -> i, PlotTheme -> "Scientific"],
  Plot[{qrsol[i][k0*z]^(-1/2)/k0, qrsol0[k0*z]^(-1/2)/k0}, {z, 0, 
    zMax}, PlotRange -> All, FrameLabel -> {"z (\[Mu]m)", "b (mm)"}, 
   PlotTheme -> "Scientific", PlotLegends -> {"sol[i]", "sol0"}]}, {i,
   3, 10}]

Note, that in the case TeV = 1000 *E^-r integrals K1,K2,K3 are diverges. From the other hand in a case TeV = 1000 *E^r we have results similar as above - see Figure 2.