# Numerical error in NIntegrate and NDsolve

Table[sbar = 100;
theta = 1;
h = 5/10;
w2 = i;
nu = 7/10;
al[t_] = sbar Sin[w2 t];
rho[t_] = 1 + nu Sin[w2 t];
r[x_, t_] = -theta (1 - x) D[rho[t], t] // Simplify;
r1[x_, t_] = -theta x (1 - x) al[t] rho[t] (1 - 2 x + h (4 x - 3)) //
Simplify;

(eqs = Subtract @@@ {D[v[x, t],
t] == -al[t] x (1 - x) D[(x + h (1 - 2 x)) v[x, t],
x] + (x (1 - x)/(2 rho[t])) D[v[x, t], x, x] + r[x, t] +
r1[x, t], v[x, 0] == 0, v[0, t] == 0, v[1, t] == 0} //
Together // Numerator // Expand // Simplify) // TableForm;

vsol = v /.
First@NDSolve[Thread[eqs == 0], v, {x, 0, 1}, {t, 0, 100 Pi/w2}(*,
MaxSteps\[Rule]10^5,MaxStepSize\[Rule]{.0002,.04},
StartingStepSize\[Rule]{.0001,.01}*), WorkingPrecision -> 30];
nint1 = NIntegrate[2 vsol[x, t], {x, 0, 1}, {t, 52 Pi/w2, 54 Pi/w2},
MaxRecursion -> 20];
Integrate[
2 (1 - x) rho[t], {x, 0, 1}, {t, 52 Pi/w2, 54 Pi/w2}]/(2 Pi/w2) +
nint1/(2 Pi/w2), {i, 10, 90, 40}]


Please suggest me how to remove the numerical errors, clearly the answer which is coming is clearly off and absurd.

• It would be helpful if you added a few more words of what you think the problem is and what you expected. Also, it might be useful for you to simplify your code. When I look at it and see a Table[...,{i,420,420,420}] I am not sure I want to spend time on this. Dec 14, 2021 at 7:26
• Sir the answer was expected to be around 2. I wrote table like this because i needed values for various i. But for i larger than 400, my answers are not correct. You can ignore this table thing because i just wanted to calculate for one particular i. Dec 14, 2021 at 8:22
• For your consideration, i have removed the table and updated the question. Dec 14, 2021 at 8:23
• You give equations with machine precision (approx. 16) and then define WorkingPrecision->30. Dec 14, 2021 at 9:30

Take the periodicity of the diffequation into account. Since t only occurs in Sin or Cos with minimal 10 t, you have peroidicity of 2 Pi/10. NDSolve for high t would be very unreliable; integrate within lower period (But not from the beginning t==0 because of disturbing effects due to initial conditions). Rationalize parameters.

sbar = 420;
theta = 1;
h = 1/2;
w2 = 10;
nu = 7/10;
al[t_] = sbar Sin[2 w2 t];
rho[t_] = 1 + nu Sin[w2 t];
r[x_, t_] = -theta (1 - x) D[rho[t], t] // Simplify;
r1[x_, t_] = -theta x (1 - x) al[t] rho[t] (1 - 2 x + h (4 x - 3)) //
Simplify;

(eqs = Subtract @@@ {D[v[x, t],
t] == -al[t] x (1 - x) D[(x + h (1 - 2 x)) v[x, t],
x] + (x (1 - x)/(2 rho[t])) D[v[x, t], x, x] + r[x, t] +
r1[x, t], v[x, 0] == 0, v[0, t] == 0, v[1, t] == 0} //
Together // Numerator // Expand // Simplify) // TableForm

vsol = v /.
First@NDSolve[Thread[eqs == 0], v, {x, 0, 1}, {t, 0, 3},
MaxSteps -> 10^5, MaxStepSize -> {.0002, .04},
StartingStepSize -> {.0001, .01}]


Periodiity in t for all x shown in graphic. Integrating from t=2 Pi/w2 to t=4 Pi/w2 should yield the same as your given range.

Manipulate[
Plot[{vsol[x, t], vsol[x, t + 2 Pi/w2], vsol[x, t + 4   Pi/w2]}, {t,
2 Pi/w2, 4 Pi/w2}, PlotPoints -> 100], {{x, 1/2}, 0, 1}]

Plot3D[vsol[x, t], {x, 0, 1}, {t, 2 Pi/w2, 4 Pi/w2}, PlotRange -> All,
PlotPoints -> 50]

nint1 = NIntegrate[2 vsol[x, t], {x, 0, 1}, {t, 2 Pi/w2, 4 Pi/w2}]

(*   -0.00674974   *)

(*   Error message: ....NIntegrate obtained -0.00674974 and 6.36785*10^-8 for the integral \
and error estimates   *)


If you are not satisfied with the estimated error of 6 10^-8, you have to use higher WorkingPrecision. But this takes very long and may crash due to not enough memory.

Integrate[
2 (1 - x) rho[t], {x, 0, 1}, {t, 42 Pi/w2, 44 Pi/w2}]/(2 Pi/w2) +
nint1/(2 Pi/w2)

(*   0.989257   *)

• Thank you so much!!!!!! but the periodicity for small t is not strong. I wanted to avoid the initial disturbances and thats why integrating over large time. I will look at it. Many thanks. Dec 15, 2021 at 6:59
• Could you also please mention that while I am integrating over x and try to plot it for time only. It is not able to integrate . Plot[{NIntegrate[vsol[x, t], {x, 0, 1}]}, {t, 2 Pi/w2, 4 Pi/w2}] Dec 15, 2021 at 7:05
• Tell NIntegrate not to try symbolicProcessing for that interpolating function vsol Plot[NIntegrate[vsol[x, t], {x, 0, 1}, Method -> {Automatic, "SymbolicProcessing" -> 0}], {t, 2 Pi/w2, 4 Pi/w2}]  Dec 15, 2021 at 7:57
• Sir, could you please look the edited one, this also contains error and give me wrong solution. Dec 15, 2021 at 11:53
• On my older version 8.0 this works well. As alternative you can try wsol[x_?NumericQ, t_?NumericQ] = vsol[x, t]; nint[t_?NumericQ] := NIntegrate[wsol[x, t], {x, 0, 1}]; Plot[nint[t], {t, 2 Pi/w2, 4 Pi/w2}]  Dec 15, 2021 at 12:10