7
$\begingroup$

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

enter image description here

Improve this question
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – xentity1x
    Jun 30, 2015 at 0:09
  • $\begingroup$ 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"? $\endgroup$
    – J. M.'s persistent exhaustion ♦
    Jun 30, 2015 at 0:27
  • $\begingroup$ 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. $\endgroup$
    – Rom38
    Feb 21, 2020 at 5:35

2 Answers 2

Reset to default
8
$\begingroup$

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

enter image description here

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

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – xentity1x
    Jun 30, 2015 at 18:00
  • $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Jun 30, 2015 at 18:07
  • $\begingroup$ @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. $\endgroup$
    – bbgodfrey
    Jul 1, 2015 at 4:25
6
$\begingroup$

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

Mathematica graphics

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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy