# Handling “ill-conditioned” system of ODE's with NDSolve

I am currently dealing with a system of coupled ODE's which I would like to solve numerically. I have already implemented the system in Mathematica using NDSolve. The system has fixed boundary values (let us say at $-1$ and $1$ the values of the functions are fixed to 0) and one parameter, $k$, which I would like to vary. The system originally includes a Dirac delta, at a certain point $u'$, so my particular implementation tries to find functions at the left and at the right of $u'$ and pastes them appropriately (including the required jump in the derivative).

For small values of $k$ the system is solved properly, however for bigger values, I am getting warnings and messages about the system being "ill-conditioned" or that the "boundary conditions may not be enough to completely specify a solution". I have tried various methods for the solver including StiffnessSwitching as well as tweaking the the working precision, the Accuracy, MaxStepSize and so on.

I actually solve the system in the interval (-.9,.9) in practice and I have noticed that the system can be solved for bigger $k$ if a start making the interval smaller and smaller, however this is not an option for my setting. I believe the instabilities are coming from the boundaries, I should mention they are originally conditions at infinity and I am using a $u=\tanh(x)$ to map $\mathbb{R}$ into $(-1,1)$. Are there any tricks or recommendations on how to handle this phenomena, or alternatives to NDSolve?

EDIT: A toy model of the problem already appears here I believe $$\delta(u-u') = \left(2u\frac{\partial}{\partial u}-(1-u^2)\frac{\partial^2}{\partial u^2}+\frac{|\vec{k}|^2}{(1-u^2)}\right)G(u)$$ this is basically a Laplacian and I would like to obtain $G(u)$ for arbitrary $k$ (at least 50, in practice).

myInf = 1;
eq1[G_Symbol, u_, k_,uP_] := (2 u*G'[u] - (1 - u^2) G''[u]) + (k^2) G[u]/(1 - u^2) == DiracDelta[u - uP];
sol[uP_, k_, tol_, vv_] := NDSolve[Evaluate[{eq1[gg, u, k, uP], gg[-myInf + tol] == 0, gg'[-myInf + tol] == vv}], gg, {u, -myInf + tol, myInf - tol}];
greenSol[uP_, k_, tol_, vv_?NumericQ] := gg[u] /. Evaluate[sol[uP, k, tol, vv]][] /. u -> myInf - tol;


Example of good result:

initD = vv /. FindRoot[Evaluate[greenSol[2/5, 2, 0.0130, vv]], {vv, 0.1},  WorkingPrecision -> MachinePrecision];
Plot[Evaluate[gg[u] /. sol[2/5, 2, 0.0130, initD]], {u, -myInf + .01,  myInf - .01}, PlotRange -> {{-myInf, myInf}, Automatic}, PlotRangePadding -> Scaled[.1]]

• I'm afraid that even with your detailed description, it is hard to give a definitive answer. Would it be possible that you provide a synthetic example that shows the same behavior and that can be posted as code here? – halirutan Jun 5 '18 at 16:00
• People here generally like users to post code as Mathematica code instead of images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. You may find this this meta Q&A helpful – Michael E2 Jun 5 '18 at 16:54
• Oh I see, I will produce a simple version shortly. – ohneVal Jun 5 '18 at 17:12
• The equation in the question can be solved symbolically. Are you attempting to obtain a good procedure for obtaining it numerically as well? – bbgodfrey Jun 6 '18 at 14:39
• Yes I know, I removed a lot more stuff from the actual system. So that we focus on the problematic part, so yes I would like a numerical procedure which I can later adapt. – ohneVal Jun 6 '18 at 16:29

## 1 Answer

Because the homogeneous ODE is singular at its endpoints, begin by determining the asymptotic solution at one of them, here u == 1.

eq = 2 u*G'[u] - (1 - u^2) G''[u] + k^2 G[u]/(1 - u^2);
AsymptoticDSolveValue[eq == 0, G, {u, 1, 1}]
(* (1 - ((k/4 + k^2/4) (-1 + u))/(1 + k/2 + 1/2 (1 + k/2) k - k^2/4))
(-1 + u)^(k/2) C
+ (1 - ((-(k/4) + k^2/4) (-1 + u))/(1 - k/2 - 1/2 (1 - k/2) k - k^2/4))
(-1 + u)^(-k/2) C *)


The solution proportional to C satisfies G == 0. However, its first k/2 - 1 derivatives also vanish at the boundary, which causes NDSolve problems for large k. It is convenient, therefore, to define a new function, g[u] == G[u]/(1 - u)^(k/2), which is better behaved there, and solve it numerically (here with uP == 1/2 and k == 50.

eqp = Simplify[(1 - u)^(-k/2) eq /. G -> Function[u, g[u] (1 - u)^(k/2)]];
stp = NDSolveValue[{eqp == 0, g[1 - 10^-5] == 1, g'[1 - 10^-5] == 0} /. k -> 50,
g, {u, 1/2, 1 - 10^-5}, WorkingPrecision -> 30];


Now, do the same at u == -1.

AsymptoticDSolveValue[eq == 0, G, {u, -1, 1}]
(* (1 + u)^(k/2) (1 - ((-(k/4) - k^2/4) (1 + u))/
(1 + k/2 + 1/2 (1 + k/2) k - k^2/4)) C
+ (1 + u)^(-k/2) (1 - ((k/4 - k^2/4) (1 + u))/
(1 - k/2 - 1/2 (1 - k/2) k - k^2/4)) C *)


which indicates that here the auxiliary function should be g[u] == G[u]/(1 + u)^(k/2).

eqm = Simplify[(1 + u)^(-k/2) eq /. G -> Function[u, g[u] (1 + u)^(k/2)]];
stm = NDSolveValue[{eqm == 0, g[-1 + 10^-5] == 1, g'[-1 + 10^-5] == 0} /. k -> 50,
g, {u, -1 + 10^-5, 1/2}, WorkingPrecision -> 30]


To create the corresponding Green's Function at uP == 1/2, match stp[u] (1 - u)^(k/2) and stm[u] (1 + u)^(k/2) and their first derivatives there.

0 == (dp stp[u] (1 - u)^(k/2) - dm stm[u] (1 + u)^(k/2)) /. {u -> 1/2, k -> 50};
-1/(1 - u^2) == D[dp stp[u] (1 - u)^(k/2) - dm stm[u] (1 + u)^(k/2), u]
/. {u -> 1/2, k -> 50};
sd = Solve[{%, %%}, {dm, dp}] // Flatten;


The Green's Function at uP == 1/2 then is

gf[u_] = Piecewise[{{dm (1 + u)^(k/2) stm[u], u < 1/2},
{dp (1 - u)^(k/2) stp[u], u >= 1/2}}] /. sd /. k -> 50;

Plot[gf[u], {u, -1 + 10^-5, 1 - 10^-5}, PlotRange -> All, AxesLabel -> {u, G}] The peak value,

gf[1/2]
(* 0.00999999999999999817834769562 *)


compares favorably with the analytical value of 1/100. To see how extreme the variation of gf is, consider

gf[1 - 10^-5]
(* 2.52543880312194881294211586*10^-123 *)


It is no wonder that NDSolve has great difficulty with the original equation. Incidentally, the symbolic Green's Function can be derived without difficulty (although not with GreenFunction, unfortunately). When plotted, it matches well the solution determined numerically here.

Addendum: What to do if AsymptoticNDSolve unavailable.

In response to a comment below, the auxiliary function also can be obtained using Series as follows. For u == 1 replace G[u] by the trial auxiliary function, f[u] (1 - u)^n, and expand eq about u == 1.

Series[(eq /. G -> Function[u, f[u] (1 - u)^n]) (1 - u)^-n, {u, 1, -1}]
(* SeriesData[u, 1, {Rational[-1, 2] k^2 f + 2 n^2 f}, -1, 0, 1] *)


This term needs to vanish, which provides the desired value for n.

SeriesCoefficient[%, -1];
Solve[% == 0, n] // Flatten
(* {n -> -(k/2), n -> k/2} *)


Select the positive solution and proceed as above.

• Hi bbgodfrey, it happens to be the case that the AsymptoticDSolveValue is only available from version 11.3 onwards. Do you know of any alternatives in previous versions? – ohneVal Jun 25 '18 at 14:51
• @ohneVal Good question. Please see the material I just added as an addendum. – bbgodfrey Jun 25 '18 at 20:33