Message Boards Message Boards

Empty result when solving two differential equation with NDSolve[ ]?

Posted 2 years ago
q[x_?NumericQ]:=(830397. - 0.375 x^(87/296))/(
     4.51833*10^8 x^(34/37) - 1.02416*10^7 x - 1. x^(383/296));
    sol1 = NDSolve[
        Rationalize[{y1'[x] == 
         y1[x] (-3 P1[x] - q[x] + G'[x]/G[x]), 
        3 P1[x] y2[x] + y2'[x] == 0.95 y1[x] q[x],
        4 P1[x] y3[x] + (y1[x] G'[x])/G[x] + 
          y3'[x] == 0,
        4 P1[x] y4[x] +y4'[x] == 0.05 y1[x] q[x],
        P1[x] == 2.372463129699238`*^-19 Sqrt[y2[x] + y4[x] + y1[x] + y3[x]],
        G'[x] == (-8.722718060807123`*^73 + 
          1.80559635200202`*^-76 G[x]^4 y3[x])/G[x]^2, 3 y1[10^-30] == 10^51 \[Pi]^2, y2[10^-30] == 0,
          y3[10^-30] == 3.3*10^56 \[Pi]^2, 
          y4[10^-30] == 0, G[10^-30] == 5.468912167778173`*^27}], {y1, 
         y2, y3, y4, G}, {x, 10^-30, 10^5}, 
        Method -> "Automatic", AccuracyGoal -> 100, 
        PrecisionGoal -> 5, WorkingPrecision -> 80]

How to solve such coupled numerical differential equation, since it always gives back some empty result instead of some Interpolating Function.

POSTED BY: John Wick
5 Replies

Maybe replace the algebraic equation by a function definition. Reduce the AccuracyGoal. I see no need for such a high setting, and it reduces the step size to zero. You could raise the PrecisionGoal if a higher quality solution is needed.

Following integrates over the whole interval:

Block[{P1},
 q[x_?NumericQ] := (830397. - 
     0.375 x^(87/296))/(4.51833*10^8 x^(34/37) - 1.02416*10^7 x - 
     1. x^(383/296));
 P1[x_] := 2.372463129699238`*^-19 Sqrt[y2[x] + y4[x] + y1[x] + y3[x]];
 sol1 = NDSolve[
   SetPrecision[{y1'[x] == y1[x] (-3 P1[x] - q[x] + G'[x]/G[x]), 
     3 P1[x] y2[x] + y2'[x] == 0.95 y1[x] q[x], 
     4 P1[x] y3[x] + (y1[x] G'[x])/G[x] + y3'[x] == 0, 
     4 P1[x] y4[x] + y4'[x] == 0.05 y1[x] q[x],
     (*P1[x]==2.372463129699238`*^-19 Sqrt[y2[x]+y4[x]+y1[x]+y3[x]],*)
     G'[x] == (-8.722718060807123`*^73 + 
         1.80559635200202`*^-76 G[x]^4 y3[x])/G[x]^2, 
     3 y1[10^-30] == 10^51 \[Pi]^2, y2[10^-30] == 0, 
     y3[10^-30] == 3.3*10^56 \[Pi]^2, y4[10^-30] == 0, 
     G[10^-30] == 5.468912167778173`*^27}, 100], {y1, y2, y3, y4, 
    G}, {x, 10^-30, 10^5}, Method -> "Automatic", AccuracyGoal -> 20, 
   PrecisionGoal -> 5, WorkingPrecision -> 80]
 ]
POSTED BY: Michael Rogers
Posted 2 years ago

Thanks a lot Michael

POSTED BY: John Wick

If you didn't suppress the error message with Quiet, you would see what it thinks the problem is: 6 equations (DEs), 7 unknowns (dependent variables).

POSTED BY: Michael Rogers
Posted 2 years ago

Now I have modified my code (except rescaling of variables, for which I am not clear how to impose this).

POSTED BY: John Wick
Posted 2 years ago

Just as a simple quick experiment, can you

scale all your 10^-30 and 10^-76 and 10^-19 up nearer to 10^-2

scale all your 10^51 and 10^73 and 10^56 and 10^27 and 10^5 down nearer to 10^2

include your definition of q[x]

possibly eliminate your use of P[x] because that doesn't seem to be the usual style

eliminate all your /G[x] denominators by multiplying

remove the Quiet and DeleteCases

and then see whether it is able to produce plausible interpolating functions?

If so then gradually and gently begin slowly backing out each of those changes a little bit at a time and confirming that it still solves?

POSTED BY: Bill Nelson
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract