# Reject diverging solution of NDSolve

I'm trying to numerically simulate a spring system with complex stiffness. In essence systems of the form

$x''(t)+ (a+ ib) x(t)=0$

For this simple example an analytic solution is easy to find. The code

Eqn = {x''[t] + (a + I*b)*x[t] == 0};
AnalyticSol = DSolve[Eqn, x[t], t]


produces the output

$\left\{x(t)\to c_1 e^{t \sqrt{-a-i b}}+c_2 e^{t \left(-\sqrt{-a-i b}\right)}\right\}$

The problem with this solution is that the first term diverges. A literature search shows that the accepted solution is to just $c_1=0$, rejecting the diverging solution. This replaces one of the initial conditions. The other initial condition is then made to be complex to account for both the initial position and velocity.

Now the actual system I'm working with is much more complex and can only be solved numerically. So I was wondering if there was any way to do something like this with NDSolve. As you can see by default it will not reject the diverging solution

a = 1.1;
b = 1.2;
T = 20;
Eqn = {x''[t] + (a + I*b)*x[t] == 0};
Inits = {x[0] == 1, x'[0] == 0};
Sys = Join[Eqn, Inits];
NumericalSol = NDSolve[Sys, x[t], {t, 0, T}];
Plot[Evaluate[Re[x[t]] /. NumericalSol], {t, 0, T}]


• Now that I think about it. That's what I'm trying to figure out. Initial conditions that generate the decaying solution. I was thinking something along the lines of x'(0)+Isqrt(a+Ib) x(0)=0. But I'm having trouble getting that to work. – xentity1x Jun 30 '15 at 0:09
• BTW: what happens if you use a different method? It is known that multistep methods (used by Mathematica as the default) sometimes have trouble keeping growing solutions off. Have you tried, say, Method -> "Extrapolation"? – J. M.'s discontentment Jun 30 '15 at 0:27
• You can use analytical {x''[t] + (a + I*b)*x[t] == 0, x[0]==1}; This gives you x[t] -> E^(-Sqrt[-a - I b] t) (1 - C[1] + E^(2 Sqrt[-a - I b] t) C[1]). Just replace the C[1] with zero from the physical meaning and fill good. – Rom38 Feb 21 at 5:35

Alternatively, this can be treated as a boundary-value problem.

a = 1.1; b = 1.2; tf = 20;
eqn = {x''[t] + (a + I*b)*x[t] == 0};
inits = {x[0] == 1, x[tf] == 0};
sys = Join[eqn, inits];
sol = NDSolve[sys, x[t], {t, 0, tf}];
Plot[Evaluate[ReIm[x[t]] /. sol], {t, 0, tf}]


MichaelE2's caveat applies here too, "It might be significantly more difficult in a complicated system".

• Do you know of any place I can read about techniques like this? One problem I notice is that the left initial condition is not satisfied in the solution as I start making the right boundary further away. – xentity1x Jun 30 '15 at 18:00
• All you would ever want to know (and more) is located here. If your answer is not accurate enough, try the "Shooting" Method, perhaps coupled with increased WorkingPrecision. – bbgodfrey Jun 30 '15 at 18:07
• @xentity1x I replaced [x20] == 0 by x[tf] == 0 in the answer for generality. With this change, a = 11/10; b = 12/10; tf = 100;, and WorkingPrecision -> 30 in NDSolve works fine. – bbgodfrey Jul 1 '15 at 4:25

If you can define an objective function that measures the size of the solution, you could optimize it. This is simple to do on the simple test case. It's a linear system, so the convergence/divergence will depend only on the ratio x'[0]/x[0]. One can optimize varying x'[t] for x[0] == 1 and test x[0] == 0 separately (best to do x[0] first, but I omit the test here). Since the system is complex, we need to consider complex-value initial conditions.

Clear[tf];
a = 1.1;
b = 1.2;
eqn = {x''[t] + (a + I*b)*x[t] == 0};
inits = {x[0] == s, x'[0] == v};
sys = Join[eqn, inits];
psol = ParametricNDSolveValue[sys, x, {t, 0, tf}, {s, v, tf}];

obj[s_?NumericQ, v_?NumericQ, tf_: 100] :=
Abs[psol[s, v, tf][psol[s, v, tf]["Domain"][[-1, -1]]]];

{min, icsol} = NMinimize[obj[1, vR + vI I, 20], {vR, vI}]
(*  {0.0000762772, {vR -> -0.513752, vI -> 1.16788}}  *)

Plot[Through[{Re, Im}[psol[1, vR + vI I /. icsol, 20][t]]], {t, 0, 20}]


It might be significantly more difficult in a complicated system, but one has to develop some measure of when one is getting close to the desired solution to be able to do this.

(Sorry about trivial the refactoring. It's best to avoid starting your own variables with capitals, since they cannot conflict with Mathematica symbols.)