Reject diverging solution of NDSolve

Question

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}]

Answer 1

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".

Answer 2

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.)