# How to solve this set of integrated PDEs

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

• This account is yours, right? It's better to merge your acount: meta.stackoverflow.com/questions/256689/… Commented Oct 15, 2021 at 10:59
• I did not know about two accounts. I submitted the request, thank you. Commented Oct 16, 2021 at 2:09
• Your problem can be solved with NDSolve[]  in a case when K1, K2, K3 can be defined as a functions of \[Nu], qr, qi in a closed form . Your first code is some mentioned case since integrals K1, K2, K3 are computed exactly . In general case we have a system of integrodifferential equations, and it is why NDSolve failed with your second code . I can recommend some FDM method to solve this problem Commented Oct 16, 2021 at 8:13
• That is great! I need the FDM method. Thank you! Commented Oct 16, 2021 at 9:04
• @AlexTrounev Could you help to recomment the FDM method for this problem? Commented Oct 19, 2021 at 13:46

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.