# Numerical integration precision issue

 theta = 1/10;
ep = 15;
w2 = 10;
c = 1/2;
rho[t_] := 1 + c Sin[  w2 t ];
h = 1/2;
S1 = ep;

sbar[t_] := ep rho[t];
S = rho[t] sbar[t];
sol = NDSolve[{D[g[x, t],
t] == (x (1 - x))/(2 rho[t]) D[g[x, t], x, x] -  (x (1 - x) ep)/
2 D[  g[x, t], x], g[x, 0] == 0  , g[1, t] == 0,
g[0, t] == theta rho[t] (1 - Exp[-1000000 t]) },
g, {t, 0, (120 Pi)/w2}, {x, 0, 1}, WorkingPrecision -> 30]
f[x_, t_] := (g /. sol[[1, 1]])[x, t];

w2/(8 Pi) NIntegrate[
f[x, t], {x, 0.00001, 1 - 0.00001}, {t, (108 Pi)/w2, 116 Pi/w2}]
theta (1. - 1/ep)


There is a numerical error in the above integral, i am sure the answer is not correct and it gives me error too, i tried using working precision and max recursion but could not resolve the error.

• Did you recognize the NDSolve messages? Oct 20, 2021 at 10:52
• I am sorry sir, i don't know how to resolve this error. Oct 20, 2021 at 11:07
• What Message did NDSolve or NIntegrate give? Oct 20, 2021 at 12:26
• @MichaelE2 NDSolveValue::eerr: Warning: scaled local spatial error estimate of 4409.8008907054445 at t = 37.6991118430775188615517205993540346103730. in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 25 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options. Oct 20, 2021 at 12:27
• Thanks @Ulrich. Then the question should be about NDSolve and NIntegrate should be removed from the question. If the NDSolve problem is fixed, then a separate question about NIntegrate may be asked if still necessary. Oct 20, 2021 at 12:36

Have a look at the boundary condition g[0, t] == theta rho[t] (1 - Exp[-1000000 t]). The part (1 - Exp[-1000000 t]) changes from zero to 1 in nearly dt=10^-6 time shift.

Try to increase dt

dt=0.01;
G = NDSolveValue[{D[g[x, t],t] == (x (1 - x))/(2 rho[t]) D[g[x, t], x, x]- (x (1 - x) ep)/2 D[g[x, t], x], g[x, 0] == 0, g[1, t] == 0,
g[0, t] == theta rho[t] (1 - Exp[-  (t/dt)])},
g, {t, 0, (120 Pi)/w2 }, {x, 0, 1} ]

Plot3D[G[t, x], {t, 0, (120 Pi)/w2}, {x, 0, 1}, PlotRange -> All,AxesLabel -> {"t", "x", "G[t,x]"}]


Solution doesn't change if you reduce dt further more!

Check result

w2/(8 Pi) NIntegrate[G[x, t], {x, 0 , 1 }, {t, (108 Pi)/w2, 116 Pi/w2}, Method -> "LocalAdaptive"]
(*0.0873137*)

• The answer is not correct ulrich. why i kept 10^-6 to match the boundary condition g(0,0) should be 0. Actually the conditions are g(1,t)=0, g(x,0)=0 and g(0,t)= theta rho[t] (1-delta(t)) Oct 22, 2021 at 9:03
• In my answer g[0,0]==0, numerical solution seems to be correct. Oct 22, 2021 at 9:54