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

Questions
, ,
1

q[x_?NumericQ]:=(830397. - 0.375 x^(87/296))/(
4.5183310^8 x^(34/37) - 1.0241610^7 x - 1. x^(383/296));
sol1 = NDSolve[
Rationalize[{y1’ ==
y1 (-3 P1 - q + G’/G),
3 P1 y2 + y2’ == 0.95 y1 q,
4 P1 y3 + (y1 G’)/G +
y3’ == 0,
4 P1 y4 +y4’ == 0.05 y1 q,
P1 == 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.

2

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?

3

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).

4

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

5

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]
 ]
6

Thanks a lot Michael