Message Boards Message Boards

GROUPS:

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

Posted 25 days ago
302 Views
|
5 Replies
|
1 Total Likes
|
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.

5 Replies
Posted 25 days 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?

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 24 days ago

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

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 24 days ago

Thanks a lot Michael

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