# nlnum problem of NDSolve

I'm having a difficulty to NDSolve a ODE due to nlnum (according to the error message). Below is the code I have. First of all are functions:

energy = Subdivide[0, 50, 2]; theta = Subdivide[0, 1, 2]; ni = Length[theta]; nj = Length[energy];
For[j = 2, j <= nj, j++,
For[i = 1, i <= ni, i++,
P[theta[[i]], energy[[j]], t] = {Px[theta[[i]], energy[[j]], t],
Py[theta[[i]], energy[[j]], t], Pz[theta[[i]], energy[[j]], t]}]];
For[j = 2, j <= nj, j++,
For[i = 1, i <= ni, i++,
Pbar[theta[[i]], energy[[j]],
t] = {Pbarx[theta[[i]], energy[[j]], t],
Pbary[theta[[i]], energy[[j]], t],
Pbarz[theta[[i]], energy[[j]], t]}]];
B={0.02,0,0.9998};
MPP = Flatten[
Table[D[P[theta[[i]], energy[[j]], t], t] ==
Cross[1/energy[[j]] B,
P[theta[[i]], energy[[j]], t]] +
Cross[Sum[(1 - theta[[ii]] theta[[i]]) (theta[[ii]] -
theta[[ii - 1]])
Sum[(energy[[jj]] P[theta[[ii]], energy[[jj]], t] -
energy[[jj]] Pbar[theta[[ii]], energy[[jj]],
t]) (energy[[jj]] - energy[[jj - 1]]), {jj, 2,
nj}], {ii, 1, ni}], P[theta[[i]], energy[[j]], t]], {i, 1,
ni}, {j, 2, nj}]];
MPV = Flatten[
Table[P[theta[[i]], energy[[j]], t] == {0, 0, 1} /. t -> 10, {i, 1,
ni}, {j, 2, nj}]];
MPbarP = Flatten[
Table[D[Pbar[theta[[i]], energy[[j]], t], t] ==
Cross[-1/energy[[j]] B,
Pbar[theta[[i]], energy[[j]], t]] +
Cross[Sum[(1 - theta[[ii]] theta[[i]]) (theta[[ii]] -
theta[[ii - 1]])
Sum[(energy[[jj]] P[theta[[ii]], energy[[jj]], t] -
energy[[jj]] Pbar[theta[[ii]], energy[[jj]],
t]) (energy[[jj]] - energy[[jj - 1]]), {jj, 2,
nj}], {ii, 1, ni}], Pbar[theta[[i]], energy[[j]], t]], {i,
1, ni}, {j, 2, nj}]];
MPbarV = Flatten[
Table[Pbar[theta[[i]], energy[[j]], t] == {0, 0, 1} /. t -> 10, {i,
1, ni}, {j, 2, nj}]];
Ms = Join[MPP, MPbarP, MPV, MPbarV];
Msss = NDSolve[Ms, Pz[theta[[2]], energy[[2]], t], {t, 10, 50},
MaxSteps -> \[Infinity]]

My code is very long and hard to read, sorry about that.

You can find that P and Pbar are vectors. I don't know is there any way better to define them. If you know please tell me.

If you have any ways to make my code easy to read(write in another way), thanks to tell me.

The ODE question I want to solve have two integrates and I have no idea how to do it and make it simple(in my thoughts), so I use two sum.

In other to make the error more clearly and the code easier to read, the code missed many coefficients, so the equation may cant't be solved numerically.

The follow one is original equation with initial value of it, the value of coefficients are included also. I have no idea how to solve it by NDSolve. If anyone have ideas how to do this?

Coefficients:

B = {0.02, 0, 0.9998};z = {0, 0, 1};zeta = 1.202;Ebare = 10;Ebaree = 15;Ebarx = 24;betae = 0.315;betaee = 0.21;betax = 0.131;u =2.13386*10^6;w = 1/0.197;\[Lambda] =(1.92893*10^11 E^(-7.24384 ArcSin[1 - 0.0000183824 t]^1.1))/t^2.4;

Equations:

D[P[theta, energy, t],
t] == {w/energy B + \[Lambda] z +
u 2/(3 zeta)
Integrate[(1 -
theta theta1) Integrate[(((
betae (betae energy1)^2)/((E^(betae energy1) +
1) Ebare) + (
betax (betax energy1)^2)/((E^(betax energy1) +
1) Ebarx)) P[theta1, energy1,
t] - ((betaee (betaee energy1)^2)/((E^(betaee energy1) +
1) Ebaree) + (
betax (betax energy1)^2)/((E^(betax energy1) +
1) Ebarx)) Pbar[theta1, energy1, t]), {energy1, 0,
50}], {theta1, Cos[ArcSin[10/t]], 1}]}\[Cross]P[theta,
energy, t];
D[Pbar[theta, energy, t],
t] == {-w/energy B + \[Lambda] z +
u 2/(3 zeta)
Integrate[(1 -
theta theta1) Integrate[(((
betae (betae energy1)^2)/((E^(betae energy1) +
1) Ebare) + (
betax (betax energy1)^2)/((E^(betax energy1) +
1) Ebarx)) P[theta1, energy1,
t] - ((betaee (betaee energy1)^2)/((E^(betaee energy1) +
1) Ebaree) + (
betax (betax energy1)^2)/((E^(betax energy1) +
1) Ebarx)) Pbar[theta1, energy1, t]), {energy1, 0,
50}], {theta1, Cos[ArcSin[10/t]], 1}]}\[Cross]Pbar[theta,
energy, t];

Initial value:

P[theta, energy, 10] == {0,
0, ((betae (betae energy)^2)/((E^(betae energy) + 1) Ebare) - (
betax (betax energy)^2)/((E^(betax energy) + 1) Ebarx))/((
betae (betae energy)^2)/((E^(betae energy) + 1) Ebare) + (
betax (betax energy)^2)/((E^(betax energy) + 1) Ebarx))};
Pbar[theta, energy, 10] == {0,
0, ((betaee (betaee energy)^2)/((E^(betaee energy) + 1) Ebaree) - (
betax (betax energy)^2)/((E^(betax energy) + 1) Ebarx))/((
betaee (betaee energy)^2)/((E^(betaee energy) + 1) Ebaree) + (
betax (betax energy)^2)/((E^(betax energy) + 1) Ebarx))}
• Definition of B is missing. – xzczd Nov 1 '18 at 11:18
• Why do functions have 3 arguments, but in the NDSolve only t? What do the original equations look like? – Alex Trounev Nov 1 '18 at 11:42
• @xzczd definition of B={0.02,0,0.9998} has been added. Thanks. – 袁子奕 Nov 1 '18 at 11:50
• @AlexTrounev the energy term and theta term has been integrate in original equation, in my function it has been summed, so the dimension of final equation is equal to dimension of energy times dimension of theta times 2 times 3: ni*nj*2*3 – 袁子奕 Nov 1 '18 at 11:54
• We already have a number of posts about vectorization of differential equations, have you checked them? BTW, you'd better put a bit more effort in improving your question, or even consider deleting this and starting a new question (it should be asked in a better way, of course), given your target is largely changed now. As to how to improve the question, for example, you may consider adding some background information, and the $\LaTeX$ form of the equation system to the post. – xzczd Nov 2 '18 at 1:12

## 1 Answer

I will answer the second question, since the first one already has a solution. We will use the method of successive approximations. My message shows only the first step. I discovered numerical instability, so I will continue after consulting with the author of the model.As a zero approximation, we use the initial data. When integrating on theta1 variable interval {theta1, Cos[ArcSin[10/t]], 1} is mapped to {x,0,1}.

B = {0.02, 0, 0.9998}; t0 = 10; tm = 50; z = {0, 0,
1}; zeta = 1.202; Ebare = 10; Ebaree = 15; Ebarx = 24; betae = \
0.315; betaee = 0.21; betax = 0.131; u = 2.13386*10^6; w =
1/0.197; \[Lambda][
t_] := (1.92893*10^11 E^(-7.24384 ArcSin[1 - 0.0000183824 t]^1.1))/
t^2.4; q[t_] := Cos[ArcSin[10/t]]; theta1 = q[t] + (1 - q[t])*x; dth =
1 - q[t];

P0 = {0, 0, ((betae (betae energy)^2)/((E^(betae energy) +
1) Ebare) - (betax (betax energy)^2)/((E^(betax energy) +
1) Ebarx))/((betae (betae energy)^2)/((E^(betae energy) +
1) Ebare) + (betax (betax energy)^2)/((E^(betax energy) +
1) Ebarx))};
Pbar0 = {0,
0, ((betaee (betaee energy)^2)/((E^(betaee energy) +
1) Ebaree) - (betax (betax energy)^2)/((E^(betax energy) +
1) Ebarx))/((betaee (betaee energy)^2)/((E^(betaee energy) \
+ 1) Ebaree) + (betax (betax energy)^2)/((E^(betax energy) +
1) Ebarx))};
(*Go back to the label*)
I1 = ((((betae (betae energy1)^2)/((E^(betae energy1) +
1) Ebare) + (betax (betax energy1)^2)/((E^(betax \
energy1) +
1) Ebarx)) P0 - ((betaee (betaee energy1)^2)/((E^(betaee \
energy1) +
1) Ebaree) + (betax (betax energy1)^2)/((E^(betax \
energy1) + 1) Ebarx)) Pbar0));
Table[f11[i] =
Interpolation[
Flatten[Table[{{t, energy},
NIntegrate[I1[[i]]*dth, {energy1, 1, 51}, {x, 0, 1},
AccuracyGoal -> 5]}, {t, t0, tm, 2}, {energy, 1, 51, 5}],
1]], {i, 1, 3}];
Table[f12[i] =
Interpolation[
Flatten[Table[{{t, energy},
NIntegrate[theta1*I1[[i]]*dth, {energy1, 1, 51}, {x, 0, 1},
AccuracyGoal -> 5]}, {t, t0, tm, 2}, {energy, 1, 51, 5}],
1]], {i, 1, 3}];
F1 = Table[
u 2/(3 zeta)*(f11[i][t, energy] - theta *f12[i][t, energy]), {i, 1,
3}];
I2 = ((((betae (betae energy1)^2)/((E^(betae energy1) +
1) Ebare) + (betax (betax energy1)^2)/((E^(betax \
energy1) +
1) Ebarx)) P0 - ((betaee (betaee energy1)^2)/((E^(betaee \
energy1) +
1) Ebaree) + (betax (betax energy1)^2)/((E^(betax \
energy1) + 1) Ebarx)) Pbar0));
Table[f21[i] =
Interpolation[
Flatten[Table[{{t, energy},
NIntegrate[I2[[i]]*dth, {energy1, 1, 51}, {x, 0, 1},
AccuracyGoal -> 5]}, {t, t0, tm, 2}, {energy, 1, 51, 5}],
1]], {i, 1, 3}];
Table[f22[i] =
Interpolation[
Flatten[Table[{{t, energy},
NIntegrate[theta1*I2[[i]]*dth, {energy1, 1, 51}, {x, 0, 1},
AccuracyGoal -> 5]}, {t, t0, tm, 2}, {energy, 1, 51, 5}],
1]], {i, 1, 3}];
F2 = Table[
u 2/(3 zeta)*(f21[i][t, energy] - theta *f22[i][t, energy]), {i, 1,
3}];
P = {P1[t], P2[t], P3[t]}; Pbar = {PB1[t], PB2[t], PB3[t]};
eq1 = D[P, t] - Cross[w/energy B + \[Lambda][t] z + F1, P]; eq2 =
D[Pbar, t] - Cross[-w/energy B + \[Lambda][t] z + F2, Pbar];
ic = {P1[t0] == 0, P2[t0] == 0,
P3[t0] == ((betae (betae energy)^2)/((E^(betae energy) +
1) Ebare) - (betax (betax energy)^2)/((E^(betax energy) +
1) Ebarx))/((betae (betae energy)^2)/((E^(betae energy) +
1) Ebare) + (betax (betax energy)^2)/((E^(betax energy) +
1) Ebarx)),
PB1[t0] == 0, PB2[t0] == 0,
PB3[t0] == ((betaee (betaee energy)^2)/((E^(betaee energy) +
1) Ebaree) - (betax (betax energy)^2)/((E^(betax energy) \
+ 1) Ebarx))/((betaee (betaee energy)^2)/((E^(betaee energy) +
1) Ebaree) + (betax (betax energy)^2)/((E^(betax energy) \
+ 1) Ebarx))};

P10 = ParametricNDSolveValue[{Table[{eq1[[j]] == 0,
eq2[[j]] == 0}, {j, 1, 3}], ic},
P1, {t, t0, tm}, {theta, energy}];

{Plot[Evaluate[Table[P10[.5, energy][t], {energy, 10, 50, 10}]], {t,
t0, tm}],
Plot[Evaluate[Table[P10[.5, energy][t], {energy, 5, 10, 1}]], {t, t0,
tm}], Plot[
Evaluate[Table[P10[.5, energy][t], {energy, 1, 4, 1}]], {t, t0,
tm}]}

Divergence occurs at the lower limit energy = 1. The author assumed the lower limit energy = 0, but in this case the integrals are not calculated due to the error 1/0. Therefore, I moved the lower limit to energy = 1, but apparently it will be necessary to increase even more. And so, we got the solution in the first step in the form of functions P0={P10[theta, energy][t],P20[theta, energy][t],P30[theta, energy][t]}, Pbar0={PB10[theta, energy][t],PB20[theta, energy][t],PB30[theta, energy][t]}. We return to the label (*Go back to the label*) and repeat.

• I fell very sorry for connect you too late. I have some questions: 1. you mentioned that there is an error 1/0, I think this can be solve by use \\$MachineEpsilon as the low limit as MMA suggested. 2.I doubt that the use of Interpolation and NIntegrate is reasonable. In my question, the integrate term "P(theta, energy, t)" is changed, but you use initial condition to "Interpolation" one, I draw the Interpolation of f11 and find that P changed to "0" very quickly so that the third term F1 is almost equal to "0" which is certainly not what I want. – 袁子奕 Nov 8 '18 at 7:31
• All in our hands. We can use any initial approximation. I just showed the method. And you can find a solution yourself. – Alex Trounev Nov 8 '18 at 8:19