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
NDSolve`u$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} = NDSolve`ProcessEquations[sys, u, {x, -L, L}];
state["NumericalFunction"]["FunctionExpression"]
Function[{x, u, NDSolve`u$2$1, NDSolve`u$2$1$2, NDSolve`u$2$1$2$3}, {NDSolve`u$2$1, NDSolve`u$2$1$2, NDSolve`u$2$1$2$3, NDSolve`u$2$1 - u NDSolve`u$2$1 - NDSolve`u$2$1$2 - NDSolve`u$2$1$2$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;
NDSolve`Iterate[newstate]
soltest = u /. NDSolve`ProcessSolutions@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]
NDSolve
isn't able to handleNIntegrate
whose integrand involves the unknown function, AFAIK. Related: mathematica.stackexchange.com/q/189698/1871 $\endgroup$Shooting
method is used, tacklingNDSolve`ProcessEquations
won't help. (Now in v12FiniteElement
can handle nonlinear BVP, too, soFiniteElement
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. ) $\endgroup$NDSolve`` ProcessEquations
? Anyway, thanks for the suggestion! $\endgroup$