# Solve ODE via NDSolveValue with NIntegrate of the Undetermined Function Being one of Terms

I'm new here. If there is anything not appropriate pls let me know.

I am currently working on a differential equation with one of which term is a integral of the variable.

$$\frac{d^2u(x)}{dx^2}=cosh(G(x))+\frac{1}{C_1}\int_{0}^{1}{u(x)sinh(G(x))dx }+C_2$$ with the boundary condition, $$u(x=0)=0$$ and $$u'(x=1)=0$$,

where $$G(x)= 2ln(\frac{1+C_3\cdot exp(-x)}{1-C_3\cdot exp(-x)})$$ and C1, C2 and C3 are system constants which could be predefined.

To show it more simply, I let $$C_1=C_2=C_3=1$$ in the code below,

G[x_] = 2 Log[(1 + Exp[-x])/(1 - Exp[-x])];

Sol =
NDSolveValue[
{
u''[x] ==
Cosh[ G[x] ] + NIntegrate[ u[x] *Sinh[ G[x] ],{x,0,1}] + 1
,u'[1] == 0., u[0] == 0.
}
, u, {x, 0, 1}, PrecisionGoal -> 10] ;


However, since u(x) in the integral have yet been solved, the numerical integration would fail with the error message show:

"The integrand {} has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,1}}"*

Some suggested by breaking the procedure of NDSolveValue in parts with the NIntegrate inserted. However, I am not sure how to do it correctly in Mathematica.

Thanks for your kindly help, I am really appreciated it!

EDIT 1

Special thanks to Tugrul Temel and Alex Trounev, who shows the singularity point at x = 0 for G[x] function. I make a adjustment follow with Alex Trounev, show below, for making the problem solvable! $$G(x)= 2ln(\frac{1+exp(-x)}{1-C_3\cdot exp(-x)})$$ where $$C_3 = 0.99$$

• You should look at the solution carefully because there may be some variables with missing numerical values. – Tugrul Temel Dec 26 '20 at 20:09
• We have singular point at x=0. Therefore we can't handle this problem with Mathematica without regularization. – Alex Trounev Dec 26 '20 at 20:13
• If we take C3<1 then there is regular solution (I have working code for this kind of problem). – Alex Trounev Dec 26 '20 at 22:57
• @TugrulTemel Yap...That is my carelessness! Actually C3 is a system constant and it would never be 1. Thanks for your suggestion! – Yuko Dec 27 '20 at 10:04
• @AlexTrounev Thanks for understanding of the problem even a ill C3 given. – Yuko Dec 27 '20 at 10:06

G[x_] = 2 Log[(1 + c3 Exp[-x])/(1 - c3 Exp[-x])]


Since the integral is a number, name it nint and solve for it later. Get anylytical solution for this reduced equation. Do indefinite integration. (I don't show the lengthy intermediate results)

usol[c1_, c2_, c3_, nint_] =
u /. First@
DSolve[{u''[x] == Cosh[G[x]] + 1/c1 nint + c2, u'[1] == 0,
u[0] == 0}, u, x]

integrand = u[x]*Sinh[G[x]] /. u -> usol[c1, c2, c3, nint] // Simplify

mint[x_] = Integrate[integrand, x]

mmii[c1_, c2_, c3_, nint_] =
Limit[mint[x], x -> 1, Direction -> 1] -
Limit[mint[x], x -> 0, Direction -> -1] // Simplify


The integral has to be equal the preassumed number nint.

nintsol[c1_, c2_, c3_] =
nint /. First@Solve[mmii[c1, c2, c3, nint] == nint, nint] // Simplify

Manipulate[
Plot[usol[c1, c2, c3, nintsol[c1, c2, c3]][x], {x, 0, 1},
PlotRange -> {{0, 1}, Automatic},
GridLines -> Automatic], {{c1, 1}, -3, 3,
Appearance -> "Labeled"}, {{c2, 1}, -3, 3,
Appearance -> "Labeled"}, {{c3, .99}, -3, 3,
Appearance -> "Labeled"}]


• It is very nice code (+1). – Alex Trounev Jan 6 at 15:52
• @AlexTrounev, You are welcome. – Akku14 Jan 6 at 17:11
• Masterpiece, still trying to understand the code, haha – Yuko Jan 7 at 15:35

In the case C3=1 we have singular solution with $$x^{-2}$$ singularity at $$x \rightarrow 0$$. Therefore we can suggest that $$C3<0$$, and in this case we have regular solution of this problem. For instance, for C3=0.99 we can get solution with using colocation method and Haar wavelets as follows

C3 = 0.99; G[x_] := 2 Log[(1 + Exp[-x])/(1 - C3 Exp[-x])];

Get["NumericalDifferentialEquationAnalysis"];
L = 1; J = 6; M = 2^J; dx = 1/2/M; xl = Table[l*dx, {l, 0, 2*M}];
xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2*M + 1}];
pp = GaussianQuadratureWeights[2*M, 0, 1]; points = pp[[All,1]]; weights = pp[[All,2]];
h1[x_] := Piecewise[{{1, 0 <= x <= 1}, {0, True}}]; p1[x_, n_] := (1/n!)*x^n;
h[x_, k_, m_] := Piecewise[{{1, Inequality[k/m, LessEqual, x, Less, (1 + 2*k)/(2*m)]},
{-1, Inequality[(1 + 2*k)/(2*m), LessEqual, x, Less, (1 + k)/m]}}, 0]
p[x_, k_, m_, n_] := Piecewise[{{0, x < k/m}, {(-(k/m) + x)^n/n!, Inequality[k/m, LessEqual, x,
Less, (1 + 2*k)/(2*m)]}, {((-(k/m) + x)^n - 2*(-((1 + 2*k)/(2*m)) + x)^n)/n!,
(1 + 2*k)/(2*m) <= x <= (1 + k)/m},
{((-(k/m) + x)^n + (-((1 + k)/m) + x)^n - 2*(-((1 + 2*k)/(2*m)) + x)^n)/n!, x > (1 + k)/m}}, 0]
f2[x_] := Sum[af[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + a0*h1[x];
f1[x_] := Sum[af[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + a0*p1[x, 1] + f10;
f0[x_] := Sum[af[i, j]*p[x, i, 2^j, 2], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + a0*p1[x, 2] + f10*x +
f00;
bc0 = {f0[0] == 0};
bc1 = {f1[L] == 0};
var = Flatten[Table[af[i, j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]];
varM = Join[{a0, f10, f00}, var]; int = Sum[weights[[i]]*(f0[z]*Sinh[G[z]] /. z -> points[[i]]),
{i, Length[points]}];
eqf[x_] := -f2[x] + Cosh[G[x]] + int + 1; eq = Flatten[Table[eqf[x] == 0, {x, xcol}]];
eqM = Join[eq, bc0, bc1];


Note, that $f2=u'',f1=u',f0=u$, eqf[x] is equation we try to solve. Since this equation is linear we use

{b, m} = CoefficientArrays[eqM, varM]

sol = -Inverse[m].b;


Finally we can plot numerical solution for 128 (solid line) and for 64 (red points) colocation points to verifier how numerical solution converges with number of points increasing. List lst64 we can prepare with the code above for J=5

lst = Table[{x,
f0[x] /. Table[varM[[s]] -> sol[[s]], {s, Length[sol]}]}, {x,
xcol}];
Show[ListLinePlot[Join[{{0, 0}}, lst], AxesLabel -> {"x", "u"},
PlotLabel -> Row[{"C3 = ", C3}]], ListPlot[lst64, PlotStyle -> Red]]
`

• This is a great solution! Thanks for your effort ! – Yuko Dec 27 '20 at 10:22
• @Yuko You are welcome! – Alex Trounev Dec 27 '20 at 11:42
• (Hi Alex, you might consider upvoting the OP's since they are new and have almost no rep yet.) – Michael E2 Dec 28 '20 at 19:00
• @MichaelE2 Thank you very much! You are right, and +1 goes at Yuko's account. – Alex Trounev Dec 28 '20 at 20:19