# Combine plots with built-in Maxima, trajectory in Sage available?

I want to combine the following plots as one output figure:

maxima('plotdf(4*y/x+x*sqrt(y),[trajectory_at,1,1],[x,0,4],[y,0,10])')
maxima('plot2d(1/4*x^4*(log(x)+2)^2,[x,0,4],[y,0,10])')


How can I achieve that? However, if there is any way to output a trajectory by Sage directly, I would be very happy. With pure Sage, I've just achieved this so far:

plot_slope_field(4*y/x+x*sqrt(y),(x,0,4),(y,0,10),headaxislength=3,headlength=3)

edit retag close merge delete

Sort by » oldest newest most voted

You can combine plots in Maxima with the command "xfun" inside plotdf e.g.

plotdf([x+y,x/2+y], [xfun,"-0.7071*x;0.7071*x"])$ In you case: plotdf(4*y/x+x*sqrt(y),[xfun,"x^4*(log(x)+2)^2/4"], [trajectory_at,1,1],[x,0,4],[y,0,10])$


you can use several functions inside the quotes, separated by semi colons.

Daniel

more

The Sage interface is the other way around; the numerical ODE solver desolve_rk4 can produce a plot:

var('x,y')
ode_rhs = 4*y/x+x*sqrt(y)
slope_field = plot_slope_field(ode_rhs, (x,0,4), (y,0,10), headaxislength=3, headlength=3, color='darkred')
exact_sol = plot(lambda x: 1/4*x^4*(log(x)+2)^2, xmin=0, xmax=4, ymin=0, ymax=10, color='blue')
integral_curve = desolve_rk4(ode_rhs, y, ivar=x, ics=[1,1], output='plot', xmax=4, ymax=10, color='yellow')
slope_field + exact_sol + integral_curve


Indeed, the plots overlap.

more

Thank you very much! I didn't know about desolve_rk4. Due to the name I assume it's about a Runge-Kutta method. Great, the two curves overlap. I've noticed that desolve itself gives a wrong output. Try this: y=function('y')(x)desolve(diff(y)==4*y/x+x*sqrt(y),y,ics=[1,1]). It gives you a minus sign instead of a plus sign in 1/4*x^4*(log(x)+2)^2 (after factorizing the output). Mathematica outputs both. Should I open a new thread about this?

( 2019-01-11 09:16:23 -0600 )edit

Yes, you can find the documentation of desolve_rk4 in Sage by typing desolve_rk4?. I found it in the online reference manual under Solving differential equations. Please open a new thread about the problem with desolve.

( 2019-01-11 10:28:16 -0600 )edit

You can add plots in Sage, so in your case:

sage: plot_slope_field(4*y/x+x*sqrt(y),(x,0,4),(y,0,10),headaxislength=3,headlength=3, color='red') + plot(1/4*x^4*(log(x)+2)^2,0,4, ymax=10, color='blue')

more

Actually, I want to check if the computed trajectory (by Maxima) coincides with my function graph. Is there a way to combine those plots to one figure via Maxima? Or alternatively, is there a way to make Sage compute trajectories like Maxima does? I need the combination of both abilities, the one by Sage (combining several plots to one) and the other by Maxima (computing/displaying trajectories without having to solve the differential equation manually).

( 2019-01-10 21:34:23 -0600 )edit

What feature of maxima is missing in Sage ? The only code from maxima that i see is the plotting, so i just explained how to move that part in Sage.

( 2019-01-11 07:23:51 -0600 )edit

Specifying trajectory_at in Maxima produces an integral curve, which you can get in Sage e.g. by desolve_rk4 (see my answer).

( 2019-01-11 08:24:10 -0600 )edit
1

"What feature of maxima is missing in Sage ?" Computing trajectories. In Maxima, after the direction field has been plotted, you can click on some point and then there will be displayed the corresponding trajectory. A great interactive feature! I haven't seen that in Sage yet.

( 2019-01-11 09:26:15 -0600 )edit

[ Not really an answer, but the damn server software refuses this as a comment... ]

I don't understand your beef with desolve's answer :

sage: y=function("y")(x)
sage: E1=diff(y(x),x)==4*y(x)/x+x*sqrt(y(x))
sage: S1=desolve(E1,y(x));S1
1/4*(2*_C + log(x))^2*x^4
sage: var("_C")
_C
True


Therefore, S1 is a set of solutions of E1 (though possibly not the set of solutions of S1, but this is another question).

BTW, both mathematica and sympy give the same answers (modulo presentation and constant's name) :

sage: mathematica("Factor[Part[DSolve[D[y[x],x]==4*y[x]/x+x*Sqrt[y[x]],y[x],x],1,1]]")
y[x] -> (x^4*(2*C[1] + Log[x])^2)/4
sage: import sympy
sage: sympy.dsolve(*[sympy.sympify(u) for u in [E1,y(x)]])_sage_()
y(x) == 1/4*(2*C1 + log(x))^2*x^4


algorithm="fricas" can't (currently) solve this ODE. Open question : quid of giac ?

more

Yes, the outputs are (seem to be) true for general solutions of this differential equation. But when it comes to initial values, at least (1,1), it gives a false output. Try this: y=function('y')(x);desolve(diff(y)==4*y/x+x*sqrt(y),y,ics=[1,1]).factor(). The output is 1/4*x^4*(log(x) - 2)^2 instead of 1/4*x^4*(log(x) + 2)^2. To get the two outputs in Mathematica, run this: DSolve[{D[y[x], x] == 4*y[x]/x + x*Sqrt[y[x]], y[1] == 1}, y[x], x].

( 2019-01-12 12:33:03 -0600 )edit

I see what you mean : your initial condition set admits both -1 and 1 as roots for _C, but your differenial equation E1 is verified only for _C=1. Mathematica gives you both solutions without warning, and Maxima picks the wrong one.

Film at 11...

BTW : sympy complains of NotImplementedError: Initial conditions produced too many solutions for constants, and throws the towel. giac seems to pick the right one.

Workaround :

• Solve the (system of) differential equation(s), getting a solution with symbolic constants (usually not declared : declare them to Sage...)
• Build a system composed of the ICS and the original differential equation(s).
• Solve this system for the symbolic constants.
( 2019-01-13 05:25:08 -0600 )edit

I consider this behavior a bug. Can you please open a ticket?

( 2019-01-13 21:39:24 -0600 )edit