2

I ran into a system of ordinary differential equations of second order and was puzzled by it for a few days. enter image description here It was originally easy to solve, but because the coupling term is second order, it cannot be solved directly by Mma. Here's the code I tried

eqs = {q[1]''[t] + q[1][t] -Sin[Pi t] (9.81 - q[1]''[t] Sin[ Pi t] - q[2]''[t] Sin[ 2*Pi t]) == 0 ,
  q[2]''[t] + q[2][t] - Sin[2*Pi t] (9.81 - q[1]''[t] Sin[ Pi t] - q[2]''[t] Sin[ 2*Pi t]) == 0}
int = {q[1][0] == 0, q[2][0] == 0, q[1]'[0] == 0, q[2]'[0] == 0};
var = {q[1], q[2]};
exp = eqs~Join~int;
NDSolve[exp, var, {t, 0, 1}]

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.

Infinity::indet: Indeterminate expression 0. ComplexInfinity ComplexInfinity encountered.

General::stop: Further output of Infinity::indet will be suppressed during this calculation.

Infinity::indet: Indeterminate expression 0. ComplexInfinity ComplexInfinity encountered.

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. I'm not sure if this system can be solved with mma, and I don't know how to solve it. I really hope you can help me.

xzczd
  • 65,995
  • 9
  • 163
  • 468
mozeq
  • 295
  • 1
  • 7

1 Answers1

2

It's because when transforming the system to the standard form required by ODE solver of NDSolve, a removable singularity is generated at $t=0$. (For more information check this post. ) There're at least 3 solutions for this problem. The easiest one is add SolveDelayed -> True:

NDSolveValue[exp, var, {t, 0, 1}, SolveDelayed -> True]

Or use a small number as starting point instead of $0$:

eps = 10^-3;
intappro = {q[1][eps] == 0, q[2][eps] == 0, q[1]'[eps] == 0, q[2]'[eps] == 0};
NDSolveValue[{eqs, intappro}, var, {t, eps, 1}]

Or use the asymptotic solution near $t=0$ as the initial condition:

asysol = AsymptoticDSolveValue[exp, var, {t, 0, 5}]

asyexpr = {q[1]@t, q[2]@t} == asysol // Thread

newic = {asyexpr, D[asyexpr, t]} /. t -> eps
sol = NDSolveValue[{eqs, newic}, var, {t, eps, 1}]

sol // ListLinePlot

enter image description here

xzczd
  • 65,995
  • 9
  • 163
  • 468
  • I almost thought this equation was impossible to solve, and I really wanted to solve it myself. It seems that I am not familiar with the NDSolve solver and the process of solving differential equations. Thank you very much for your help and make me saw this problem clearly. – mozeq Nov 21 '19 at 08:37