# Nonlinear boundary value problem of ODE involving principal value of integral

I have asked a question on a nonlinear eigenvalue problem (EVP). And I have worked on these for a week but I cannot solve it. I think I should first try to solve the related ODE boundary value problem on $$x\in[-L,L]$$:

$$-u^\prime+uu^\prime+u^{\prime\prime}+u^{\prime\prime\prime}+u^{\prime\prime\prime\prime}+\frac{1}{2L}PV\int_{-L}^L u^{\prime\prime\prime}(s)\cot\left(\frac{\pi(x-s)}{2L}\right)\mathrm{d}s=0,$$

where derivatives are denoted by primes, PV represents the principal value of the integral due to the singularity at $$x=s$$. It should be solved with periodic boundary condition(s): $$u(\pm L)=0$$ (and $$du/dx\vert_{x=\pm L}=0$$ if necessary) or $$u[-L]=u[L]$$. Here, the highest derivative is 4th-order, I think four BCs are needed but not very sure which BCs should be used.

This code does not work. I believe the reason is due to the NIntegrate. Can anyone give me some suggestion? Thanks in advance!

L = 20;
sys = {-u'[x] + u[x]*u'[x] + u''[x] + u'''[x] + u''''[x] +
1/(2*L)*NIntegrate[u'''[s]*Cot[(π*(x - s))/(2*L)], {s, -L, x, L},
Method -> "PrincipalValue"] == 0,
u[-L] == 0, u[L] == 0, u'[-L] == 0, u'[L] == 0};

sol = First[u /. NDSolve[sys, u, {x, -L, L}]]

It gives several warnings like this:

NIntegrate::inumr: The integrand NDSolveu$1$2$3$2792 Cot[1/40 π (-20.-s)] has evaluated to non-numerical values for all sampling points in the region with boundaries {{-20.,20}}. >>

Update: try the method provided for this question by @xzczd

L = 20;
fakefunc[x_?NumericQ] = 0;
sys = {u[x]*u'[x] + u''[x] + u''''[x] + 1/(2*L)*fakefunc[x] + u'''[x] == u'[x], u[-L] == 0, u[L] == 0, u'[-L] == 0, u'[L] == 0};

{state} = NDSolveProcessEquations[sys, u, {x, -L, L}];

state["NumericalFunction"]["FunctionExpression"]
Function[{x, u, NDSolveu$2$1, NDSolveu$$2$$1$$2, NDSolveu$$2$$1$$2$$3}, {NDSolveu21, NDSolveu212, NDSolveu2123, NDSolveu21 - u NDSolveu21 - NDSolveu212 - NDSolveu212$$3 - fakefunc[x]/40}]
Clear@int;
int[y_List, x_] := With[{func = ListInterpolation[D[y, {x, 3}], {{-L, L}}]},
NIntegrate[func@s*Cot[π*(x - s)/(2*L)], {s, -L, x, L},
Method -> {"PrincipalValue", SymbolicProcessing -> 0}]];

rule = HoldPattern@MapThread[Function[x, fakefunc[x]], {x}, 1] :> (int[u, #] & /@ x)

newstate = state /. rule;
NDSolveIterate[newstate]
soltest = u /. NDSolveProcessSolutions@newstate

which do not give any warning. But the plot gives null

Plot[Evaluate[soltest[x]], {x, -L, L}, PlotRange -> All]

Update(May 26, 2019): a probably failed case

In the following, @xzczd's answer is tested with another analytic solution $$u(x)=sech(x)$$, which has the same b.c.s as the original question.

L = 30; domain = {-L, L};
points = 61; (*could be reduced to 41*)
difforder = 4; grid = Array[# &, points, domain];
(* find pdetoae in the link above. *)

ptoafunc = pdetoae[u[x], grid, difforder];
Clear@nint
nint[u_List, x_List] := nint[u, #] & /@ x;
nint[u : {__?NumericQ}, x_] := With[{func = ListInterpolation[u, {domain}]},
NIntegrate[func[s]*Cot[(\[Pi] (x - s))/(2 L)], {s, -L, x, L},
Method -> {PrincipalValue, SymbolicProcessing -> 0}, AccuracyGoal -> 8]]

rhslst = Table[-u'[x] + u[x]*u'[x] + u''[x] + u'''[x] + u''''[x] +
1/(2*L)*NIntegrate[u'''[s]*Cot[(\[Pi]*(x - s))/(2*L)], {s, -L, x, L},
Method -> {PrincipalValue, SymbolicProcessing -> 0},
AccuracyGoal -> 8] /. u -> Function[x, Sech[x]] // Evaluate, {x, grid}];
rhsfunc = ListInterpolation[rhslst, grid];

{neweq, newbc} = {-u'[x] + u[x] u'[x] + u''[x] + u'''[x] + u''''[x] +
1/(2 L) nint[u'''[x], x] == rhsfunc@x,
With[{lhs = {u[-L] == 0, u[L] == 0, u'[-L] == 0, u'[L] == 0}[[All, 1]]},
lhs == (lhs /. u -> Function[x, Sech[x]]) // Thread]};

del = #[[3 ;; -3]] &;
ae = ptoafunc@neweq // del;
aebc = ptoafunc@newbc;
guess[x_] = 0;
solrule = FindRoot[{ae, aebc},
Table[{u@x, guess@x}, {x, grid}]]; // AbsoluteTiming
(* {2881.481811, Null} *)

solfunc = ListInterpolation[solrule[[All, -1]], grid];
Plot[solfunc@x, {x, -L, L}, PlotStyle -> Thick, PlotRange -> All]

Which is totally different the desired solution:

Plot[Sech[x], {x, -L, L}, PlotRange -> {All, {-1, 1}}, Frame -> True, Axes -> False]

• Currently NDSolve isn't able to handle NIntegrate whose integrand involves the unknown function, AFAIK. Related: mathematica.stackexchange.com/q/189698/1871 – xzczd May 24 '19 at 17:32
• Notice the method in the previous link won't fully apply to your problem, because the integration exists in the spatial dimension of a PDE system in that question, while yours is a BVP of ODE i.e. the integration exists in time dimension. As long as Shooting method is used, tackling NDSolveProcessEquations won't help. (Now in v12 FiniteElement can handle nonlinear BVP, too, so FiniteElement might help. But I haven't looked into this new feature yet. ) We probably need to build our own solver, which won't be that easy, I'm afraid. (Nonlinear BVP is already troublesome. ) – xzczd May 25 '19 at 4:38
• @xzczd oh, formidable nonlinear BVP. Previously, i just think of the difficulty of nonlinear PDE definitely larger than nonlinear BVP. If it becomes a PDE somehow, could it be handled by NDSolve` ProcessEquations? Anyway, thanks for the suggestion! – jsxs May 25 '19 at 5:25
• Yes, if the $x$ becomes the spatial dimension of a PDE, the method in the link above will be fully applicable. (Whether it'll work well or not is another story, though. My limited personal experience is, if the nonlinear BVP of ODE is hard, the corresponding PDE won't be easy to solve, either. ) – xzczd May 25 '19 at 6:00

As already mentioned in the comment, currently NDSolve isn't able to handle NIntegrate whose integrand involves the unknown function, AFAIK. The most realistic choice in my view is to leave NDSolve alone and build our own solver. The following is my solution. Though it's rather slow, it does work.

It's easy to notice the original BVP owns a trivial solution

$$u=0$$

It might owns non-trivial solution, but to make this answer more interesting and the result easier to check, I'll slightly modify the original equation and b.c.s to

$$-u'+uu'+u''+u'''+u''''+\frac{1}{2L}PV\int_{-L}^L u'''(s)\cot\left(\frac{\pi(x-s)}{2L}\right)\mathrm{d}s=\color{red}{f(x)}$$

$$u(-L)=u(L)=0$$ $$u'(-L)=u'(L)=\color{red}{-\frac{\pi}{L}}$$

where $$f(x)$$ is a properly chosen function so that the BVP owns the non-trivial solution

$$u=\sin(\frac{\pi x}{L})$$

The idea of this solution is simple: discretize the ODE and b.c.s based on finite difference method (FDM), and solve the resulting nonlinear algebraic system with FindRoot. I'll use pdetoae for the generation of difference equations.

L = 20; domain = {-L, L};
points = 50; difforder = 4; grid = Array[# &, points, domain];
(* Definition of pdetoae isn't included in this answer,
please find it in the link above. *)
ptoafunc = pdetoae[u[x], grid, difforder];
Clear@nint
nint[u_List, x_List] := nint[u, #] & /@ x;
nint[u : {__?NumericQ}, x_] :=
With[{func = ListInterpolation[u, {domain}]},
NIntegrate[func[s]*Cot[(π (x - s))/(2 L)], {s, -L, x, L},
Method -> {PrincipalValue, SymbolicProcessing -> 0}, AccuracyGoal -> 8]]

rhslst = Table[-u'[x] + u[x]*u'[x] + u''[x] + u'''[x] + u''''[x] +
1/(2*L)*NIntegrate[u'''[s]*Cot[(π*(x - s))/(2*L)], {s, -L, x, L},
Method -> {PrincipalValue, SymbolicProcessing -> 0}, AccuracyGoal -> 8] /.
u -> Function[x, Sin[Pi x/L]] // Evaluate, {x, grid}];
rhsfunc = ListInterpolation[rhslst, grid];

{neweq, newbc} = {-u'[x] + u[x] u'[x] + u''[x] + u'''[x] + u''''[x] +
1/(2 L) nint[u'''[x], x] == rhsfunc@x,
With[{lhs = {u[-L] == 0, u[L] == 0, u'[-L] == 0, u'[L] == 0}[[All, 1]]},
lhs == (lhs /. u -> Function[x, Sin[Pi x/L]]) // Thread]};
del = #[[3 ;; -3]] &;
ae = ptoafunc@neweq // del;
aebc = ptoafunc@newbc;
guess[x_] = 0;
solrule = FindRoot[{ae, aebc}, Table[{u@x, guess@x}, {x, grid}]]; // AbsoluteTiming
(* {1418.56, Null} *)
solfunc = ListInterpolation[solrule[[All, -1]], grid];
Plot[solfunc@x, {x, -L, L}, PlotRange -> All]

The right half of the solution slightly veers off the analytic solution, but I think it's not too bad. Tackling points, difforder, etc. may improve the result, but I'd like to stop here.

# Remark

1. The nint i.e. calculation for the principal value is probably the bottleneck of my solution. I've simply lived with the combination of ListInterpolation and NIntegrate, because I know nothing about numeric calculation of principal value.

2. It's lucky that FindRoot (almost) manages to find the desired result under such a casual initial guess guess[x_] = 0 in this case, but do remember looking for proper initial guess for nonlinear BVP can be painful. Just search in this site, we already have dozens of examples.

3. The original problem you want to solve i.e. nonlinear eigenvalue problem is even harder, so, grit your teeth. Another example I can remember at the moment: NDSolve eigenvalue problem of bound state

As to the test for $$u(x)=\text{sech}(x)$$ in the question, the grid chosen by OP turns out to be too coarse. After adjusting the parameters to L = 7; points = 50;, a pretty good result is obtained after 973.142 seconds. (Actually L = 7; points = 25 is good enough, and the corresponding timing is 127.755. )

• just back to my computer. I appreciate your effort very much! Please give me a few days to understand your idea, test on my real problem, and wait for potential improvement or answer from someone. I will accept one later. – jsxs May 25 '19 at 14:57
• @jsxs Feel free to wait for better answers :) . – xzczd May 25 '19 at 15:16
• Great! your idea is basically for the validation of a non-trivial solution of an inhomogeneous ODE via an "inverse problem". In other words, specify an analytic special solution then determine the inhomogeneous function via rhsfunc = ListInterpolation[...], then feed rhsfunc back to the ODE to (inversely) search the special solution iteratively with FindRoot. Am I right? – jsxs May 26 '19 at 4:49
• @jsxs Yes, exactly. – xzczd May 26 '19 at 5:03
• when solving the original homogeneous BVP with eq =... == 0; bc={u[-L] == 0, u[L] == 0, u'[-L] == 0, u'[L] == 0}; only the trival solution is found. Is it due to the initial guess? Another detailed issue, we have 4 b.c.s, thus 4 discrete b.c.s should be removed from both end, but odebc = ptoafunc@newbc not do this. Is it proper? Thank you! – jsxs May 26 '19 at 5:31

This is not a perfect answer, but I will post it for reference.

I imposed boundary conditions according to wikipedia, and tried to solve it while changing the velocity at $$x=-L$$.

First of all, since it takes time to integrate Cot [x], this function is approximated by splines.

L = 20;
cot = Interpolation[
If[Element[#2, Reals], {#1, #2}, Nothing] & @@ # & /@
Table[{x, Cot[x]}, {x, -10, 10, 0.01}]];

int[func_, x_] :=
Fold[(#1 + (1/(2*L))*func'''[#2]*cot[(\[Pi]*(x - #2))/(2*L)]) &, 0,
Range[-L, L, 0.1]];

Next, sampling points uniformly and approximate the solution u[x]
(function approximation by hermite polynomial)

u = Interpolation[Table[{i, v[i]}, {i, -L, L, 1}],
InterpolationOrder -> 4];

In this way, we convert the PDE into the optimization problem solving for the parameter v[i],which minimize $$L^2$$ norm.

Next we define the loss function, here,L^2 norm + penalty for boundary conditions.

loss[cf_?VectorQ][var_][speed_] :=
Block[{L = 20, penalty = 10000,
rule = Thread[Rule[var, cf]]},
ff = u /. rule;(*insert params*);
(*b.c.s*);
cnds = {ff[-L], ff[L] - ff[-L], ff'[-L] - speed, ff'[L] - ff'[-L],
ff''[L] - ff''[-L], ff'''[L] - ff'''[-L], ff''''[L] - ff''''[-L]};

(*L2 norm and penalty*)
Fold[(#1 +
Power[-ff'[#2] + ff[#2]*ff'[#2] + ff''[#2] + ff'''[#2] +
ff''''[#2] + 1/(2*L)*int[ff, #2], 2]) &, 0,
Range[-L, L, 1]] + Fold[#1 + penalty*Power[#2, 2] &, 0, cnds]]

Then use NMinimize to find the parameter v[i] that minimizes the loss function.

var = Table[v[i], {i, -20, 20, 1}];
sol = Last@
NMinimize[loss[var][var][speed], var, Method -> "SimulatedAnnealing"];
u[x] /. sol, {speed, Range[-1, 1, 0.1]}];

Animate[
Dynamic@Plot[answers[[m]], {x, -20, 20},
PlotLabel -> "u'[-L]=" <> ToString[((m/10) - 1.1)]], {m, 1, 21, 1}]

This is the solution we found.

From here onwards I need to understand the problem itself, so I'll leave it. I think xzdczd's method is the best.

3) For loss function,
I used L^2 norm $$||f||_2^2$$ for the part of differential equation
and penalty terms for the part of boundary condition.

In normal setting,
I have to increase the coefficient of that term,100000,every iteration,like 10x times per iteration.
but this problem needs a lot time to calcutation for one time,
so instead of update the number,I just set the big number.

if we have a true solution,the term converge to 0,so no matter how large the coffieicient,it is no problem.
that is my opinion.
(yes,my heurstics)

• My personal experience is, NMinimze, FindMinimum etc. usually won't work well on nonlinear equation solving, FindRoot seems to be the only reasonable choice in most cases. (As we can see, NMinimize fails to find the trivial solution $u=0$ in this case. ) – xzczd May 25 '19 at 9:57
• @Xminer I appreciate your effort very much! Please give me a few days to understand your idea, test on my real problem, and wait for potential improvement or answer from someone. I will accept one later. – jsxs May 25 '19 at 14:57
• @xzzczd you are correct. However, this time NMinimize itself is not problem, but it seems I have to optimize where the knots/nodes are placed. – Xminer May 25 '19 at 17:07
• @jsxs the idea is pseudo-spectral method.I wanted to use it because I'm studying it now xp. I think xzczd's answer is better. – Xminer May 25 '19 at 17:15
• @Xminer 1. I know something about pseudo-spectral, which approximates the solution on a periodic interval as a Fourier series $u(x)=\sum a_n e^{inx}$ where $a_n$ is the discrete Fourier coefficient, then the series will be truncated and the (space) derivatives can be evaluated by FFT. But I see no FFT in the code... which part is the core of pseudo-spectral? – jsxs May 26 '19 at 8:02