0
$\begingroup$

I am trying to modify the existing code developed by Michael E2 in this question here. His solution was for one differential equation. I like his code because it has ability to solve nonlinear differential equation symbolically. I would like to extend it to be able to solve four coupled differential equation of second order. I am running the following test codes and my attempt to modify the existing code:

ClearAll[x, t, a, b, c, xx];
(*******       testing for simple pendulum equation   ***************)
(*******           Series for one SONDE                 *************)
Clear[seriesDSolve];
seriesDSolve[ode_,y_,iter:{x_,x0_,n_},ics_:    {}]:=Module[{dorder,coeff},dorder=Max[0,Cases[ode,Derivative[m_][y][x]:>m,Infinity]];
(coeff=Fold[Flatten[{#1,Solve[#2==0/.#1,Union@Cases[#2/.#1,Derivative[m_][y][x0]/;m>=dorder,Infinity]]}]&,ics,Table[SeriesCoefficient[ode/.Equal->Subtract,{x,x0,k}],{k,0,Max[n-   dorder,0]}]];
Series[y[x],iter]/.coeff)/;dorder>0]

seriesDSolve[
y''[x] + a*Sin[y[x] == 0], y, {x, 0, 8}, {y[0] -> c, y'[0] -> 0}

and the output is

$c-\frac{1}{2} x^2 (a \sin (c))+\frac{1}{24} a^2 x^4 \sin (c) \cos (c)+\frac{1}{720} x^6 \left(3 a^3 \sin ^3(c)-a^3 \sin (c) \cos ^2(c)\right)+\frac{a x^8 \left(a^3 \sin (c) \cos ^3(c)-33 a^3 \sin ^3(c) \cos (c)\right)}{40320}+O\left(x^9\right)$

which is correct.

Now a test code for 4 coupled differential equations (it should also handle nonlinear ones as well).

 ClearAll[k,a];
 solution = DSolve[{x''[t] + a*x[t] + k*v[t] == 0,
                                y''[t] + a*v[t] + k*x[t] == 0, 
  u''[t] + a*v[t] + k*x[t] == 0,
  v''[t] + a*v[t] + k*x[t] == 0,
  x[0] == 1, y[0] == 1, u[0] == 1, v[0] == 1,
  x'[0] == 0, y'[0] == 0, u'[0] == 0, v'[0] == 0},
  {x, y, u, v}, {t, 0, 100}]

and the answer is :

$\left\{\left\{u\to \left(\{t\} {f4a1}\frac{1}{2} e^{t \left(-\sqrt{-a-k}\right)} \left(e^{2 t \sqrt{-a-k}}+1\right)\right),v\to \left(\{t\} {f4a1}\frac{1}{2} e^{t \left(-\sqrt{-a-k}\right)} \left(e^{2 t \sqrt{-a-k}}+1\right)\right),x\to \left(\{t\} {f4a1}\frac{1}{2} e^{t \left(-\sqrt{-a-k}\right)} \left(e^{2 t \sqrt{-a-k}}+1\right)\right),y\to \left(\{t\} {f4a1}\frac{1}{2} e^{t \left(-\sqrt{-a-k}\right)} \left(e^{2 t \sqrt{-a-k}}+1\right)\right)\right\}\right\}$

one can simplify to get trigonometric expression.

 p = FullSimplify[{x[t], y[t], u[t], v[t]} /. solution[[1]]] 

Then plotting yields,

 a = 0.2; k = 0.5;
 GraphicsGrid[{
{Plot[{x[t] /. solution[[1]]}, {t, 0, 10}], 
 Plot[{y[t] /. solution[[1]]}, {t, 0, 
 10}]}, {Plot[{u[t] /. solution[[1]]}, {t, 0, 10}], 
 Plot[{v[t] /. solution[[1]]}, {t, 0, 10}]}}]

enter image description here

Then working on series for 4 coupled SODE's,

Clear[seriesDSolve];
seriesDSolve[ode_, f_, iter : {t_, t0_, n_}, ics_: {}] := 
Module[{dorder, coeff}, 
dorder = Max[0, Cases[ode, Derivative[m_][f][t] :> m, Infinity]];
(coeff = 
 Fold[Flatten[{#1, 
     Solve[#2 == 0 /. #1, 
      Union@Cases[#2 /. #1, Derivative[m_][f][t0] /; m >= dorder, 
        Infinity]]}] &, ics, 
  Table[SeriesCoefficient[
    ode /. Equal -> Subtract, {t, t0, k}], {k, 0, 
    Max[n - dorder, 0]}]];
  Series[f[t], iter] /. coeff) /; dorder > 0]

Unfortunately, the modified code does not work. My list approach is not correct. Any suggestions on what I am missing? (MMA v11)

$\endgroup$
3
  • 2
    $\begingroup$ Why not use AsymptoticDSolveValue? $\endgroup$
    – Carl Woll
    Jan 30 at 1:59
  • $\begingroup$ @CarlWoll, AsymptoticDSolveValue is available from v12. $\endgroup$
    – Aschoolar
    Jan 30 at 12:08
  • $\begingroup$ Actually, it is available starting in 11.3. $\endgroup$
    – Carl Woll
    Jan 30 at 23:56

Your Answer

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

Browse other questions tagged or ask your own question.