1
$\begingroup$

The equation I want to solve numerically with WolframCloud is:

z[t_] := Exp[2*Sqrt[L/3]*t]
Eq1 = -40/9*L-20/3*y''[t]/y[t]+2/3*y'[t]^2/y[t]^2-6*z'[t]*y'[t]/(z[t]*y[t])+11/6*z'[t]^2/z[t]^2-360/7*1/y[t]^2+97/21*1/(y[t]^2*z[t])==0

I want to solve it for y[t] only, and then verify that it gives y[t] such that if I want to solve this equation for z[t], it returns the correct exponential function. To do so, I use the following code:

Clear[F,R,y]
L = 0.038;
sol0 = NDSolve[{Eq1, y[0]==0.1, y'[0]==5.5},y[t],{t, 0, 1}];
R[t_]:=Evaluate[y[t]/.sol0]
Eq2 =-40/9*L-20/3*R''[t]/R[t]+2/3*R'[t]^2/R[t]^2-6*F'[t]*R'[t]/(F[t]*R[t])+11/6*F'[t]^2/F[t]^2-360/7*1/R[t]^2+97/21*1/(R[t]^2*F[t])==0;
sol = NDSolve[{Eq2,F[0]==1},F[t],{t,0,10},MaxSteps->10000][[1]];
Plot[{Log[Evaluate[F[t]/.sol]], 2*Sqrt[L/3]*t},{t,0,1}]

As you see, I plot only the log of the function so it should give two superimposed lines. But it gives me this: enter image description here

We can see that until $\sim 0.6$ it gives the correct behavior, but I can't explain what is going on with my code to give a wrong answer... Have I done something wrong in my code?

EDIT: After having played with the parameters, it seems that the point where F[t] deviates from the line $2 \sqrt{\frac{L}{3}} t$ is exactly the middle of the domain of definition of R[t]. In fact, R[t] has roughly the shape of $-at^2+b$, and so its domain of definition has a maximum right at its half. I don't know what it means but if anyone has an idea...

Improve this question
$\endgroup$
7
  • $\begingroup$ Well, one obvious issue is that you solve first ode for up to $t=1$. Then you used the solution in the second ode which you now solved for up to $t=10$. So it will extrapolate the first solution, which is not good. Try to solve first ode for $t=10$. But it is stiff ode. So you can try reducing AccuracyGoal and add Method -> {"StiffnessSwitching"}, to help., Your second ode is also stiff. May be more special options are needed. I get in V 13.0.1 NDSolve::ndsz: At t == 1.3637473398450544, step size is effectively zero; singularity or stiff system suspected. $\endgroup$
    – Nasser
    Commented Feb 23, 2022 at 8:04
  • $\begingroup$ @Nasser Thank you for your comment! I've tried to change to $t=10$ in the first equation and to add the Method, but it gives the same plot at the end... $\endgroup$
    – Jeanbaptiste Roux
    Commented Feb 23, 2022 at 8:09
  • $\begingroup$ You can simplify your code a little by writing R = NDSolveValue[{Eq1, y[0] == 0.1, y'[0] == 5.5}, y, {t, 0, 1}] and now you do not need to do R[t_]:=Evaluate[y[t]/.sol0] $\endgroup$
    – Nasser
    Commented Feb 23, 2022 at 8:21
  • $\begingroup$ @Nasser I tried to do so, but it gives me errors... $\endgroup$
    – Jeanbaptiste Roux
    Commented Feb 23, 2022 at 8:57
  • $\begingroup$ @Nasser I cannot open the notebook because I didn't subscribe to a plan in WolframCloud, but did you get the right straight lines in the plot? $\endgroup$
    – Jeanbaptiste Roux
    Commented Feb 23, 2022 at 9:19

1 Answer 1

Reset to default
4
$\begingroup$ 
sol0 = NDSolve[{Eq1, y[0] == 0.1, y'[0] == 5.5}, y[t], {t, 0, 1}, InterpolationOrder -> All];

solves the problem. With the rest of the code as in the question, the final plot is

enter image description here

The difficulty with the original code is that R[t] is a third order interpolation function. Hence, R''[t] is only first order, which apparently is not accurate enough to recover z[t]. The option InterpolationOrder -> All leads to more accurate interpolation of higher derivatives.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.