Complicated System of Second Order Differential Equations

Question

I am trying to solve a coupled system of second order differential equations using either DSolve[] or NDSolve[]. I believe Mathematica may not be powerful enough to find the solution because when I use the same syntax and simpler equations, it solves them correctly but when I run the full code below, it gives me my input back.

Here is my code:

SetAttributes[epsilon, {Constant, Protected}]
SetAttributes[alpha, {Constant, Protected}]
SetAttributes[beta, {Constant, Protected}]

eqns = {phi''[t] == 
            -(3 (a[t]^2 a'[t] phi'[t] + 
            epsilon a[t]^2 a'[t] phi'[t] + 
            epsilon a[t]^2 Abs[alpha]^2 a'[t] phi'[t] + 
            epsilon Abs[alpha]^2 b[t]^2 b'[t] phi'[t]))/(a[t]^3 + epsilon a[t]^3 + 
            epsilon a[t]^3 Abs[alpha]^2 + epsilon Abs[alpha]^2 b[t]^3), 
        chi''[t] == 
            -(3 (epsilon a[t]^2 Abs[beta]^2 a'[t] chi'[t] + 
            b[t]^2 b'[t] chi'[t] + epsilon b[t]^2 b'[t] chi'[t] + 
            epsilon Abs[beta]^2 b[t]^2 b'[t] chi'[t]))/(epsilon a[t]^3 Abs[beta]^2 + 
            b[t]^3 + epsilon b[t]^3 + epsilon Abs[beta]^2 b[t]^3), 
        a''[t] == 
            1/(4 a[t]) (-2a'[t]^2 + epsilon a[t]^2 Abs[beta]^2 chi'[t]^2 + 
            a[t]^2 phi'[t]^2 + epsilon a[t]^2 phi'[t]^2 + 
            epsilon a[t]^2 Abs[alpha]^2 phi'[t]^2), 
        b''[t] == 
            1/(4 b[t]) (-2 b'[t]^2 + b[t]^2 chi'[t]^2 + 
            epsilon b[t]^2 chi'[t]^2 + epsilon Abs[beta]^2 b[t]^2 chi'[t]^2 + 
            epsilon Abs[alpha]^2 b[t]^2 phi'[t]^2)};

sol = DSolve[eqns, {phi[t], chi[t], a[t], b[t]}, t]

I've also tried putting in initial conditions and get the same result. Is it possible Mathematica just can't handle a system this complicated? Is there anything I can do to make it work? Thanks!

Answer 1

Here's a start with NDSolve. As Michael stated, need to assign values to constants and set up some initial conditions. You can change them as you like although I had to adjust them a bit before I was able to successfully integrate from t=0 to t=2 without encountering a stiff system suggesting a singularity was reached such as division by zero:

epsilon = 1;
alpha = 1.5;
beta = -2.25;

initVals = {phi[0] == 1/10, phi'[0] == 1/4, chi[0] == 1, chi'[0] == 1,
    a[0] == 1/4, a'[0] == 1/4, b[0] == -1, b'[0] == -1};
eqns = {phi''[
     t] == -(3 (a[t]^2 a'[t] phi'[t] + epsilon a[t]^2 a'[t] phi'[t] + 
          epsilon a[t]^2 Abs[alpha]^2 a'[t] phi'[t] + 
          epsilon Abs[alpha]^2 b[t]^2 b'[t] phi'[t]))/(a[t]^3 + 
       epsilon a[t]^3 + epsilon a[t]^3 Abs[alpha]^2 + 
       epsilon Abs[alpha]^2 b[t]^3), 
   chi''[t] == -(3 (epsilon a[t]^2 Abs[beta]^2 a'[t] chi'[t] + 
          b[t]^2 b'[t] chi'[t] + epsilon b[t]^2 b'[t] chi'[t] + 
          epsilon Abs[beta]^2 b[t]^2 b'[t] chi'[t]))/(epsilon a[
          t]^3 Abs[beta]^2 + b[t]^3 + epsilon b[t]^3 + 
       epsilon Abs[beta]^2 b[t]^3), 
   a''[t] == 
    1/(4 a[t]) (-2 a'[t]^2 + epsilon a[t]^2 Abs[beta]^2 chi'[t]^2 + 
       a[t]^2 phi'[t]^2 + epsilon a[t]^2 phi'[t]^2 + 
       epsilon a[t]^2 Abs[alpha]^2 phi'[t]^2), 
   b''[t] == 
    1/(4 b[t]) (-2 b'[t]^2 + b[t]^2 chi'[t]^2 + 
       epsilon b[t]^2 chi'[t]^2 + 
       epsilon Abs[beta]^2 b[t]^2 chi'[t]^2 + 
       epsilon Abs[alpha]^2 b[t]^2 phi'[t]^2)};
{f1, f2, f3, f4} = 
 NDSolveValue[Join[{eqns, initVals}], {phi, chi, a, b}, {t, 0, 2}]
Plot[{f1[t], f2[t], f3[t], f4[t]}, {t, 0, 2}]