# Plot solution for y' + 2xy = 1

f = desolve(diff(y,x) + 2xy - 1, y, ics=[0,0]); f plot(f) # error message ... unable to simplify to float approximation.

Tried plotting real_part of solution: -1/2Isqrt(pi)e^(-x^2)erf(I*x), but get same error message.

Tried using list_plot, but error about symbolic expression. Haven't been able to get implicit_plot or sol.simplify_full to work.

Thanks for any hints.

edit retag close merge delete

Sort by » oldest newest most voted

Your solution is -1/2*I*sqrt(pi)*e^(-x^2)*erf(I*x) which Sage does produce correctly. But, Sage appears to be having trouble with erf(I*x). This should be a pure imaginary result, but Sage is not giving that for $x>1.42$.

Wolfram alpha plots your function nicely and also handles erf(I*x) well. (See this link.)

Anyone have ideas on this?

more

Thank you very much for finding this!

( 2012-07-02 03:59:24 -0500 )edit

( 2012-07-02 04:42:08 -0500 )edit

You can get the plot using a numerical approximation using desolve_rk4.

y=function('y',x)
f(x,y)=1-x*y
soln=desolve_rk4(f(x,y),y,ics=[0,0],ivar=x,end_points=20)
line(soln)


Edit:

You can also plot a slope field with the solution as follows.

y=function('y',x)
f(x,y)=1-x*y
soln=desolve_rk4(f(x,y),y,ics=[0,0],ivar=x,end_points=20)
p=line(soln)
p+=plot_slope_field(f(x,y),(x,0,20),(y,0,0.8))
show(p)


Or, if you prefer a slope field with arrows.

y=function('y',x)
f(x,y)=1-x*y
soln=desolve_rk4(f(x,y),y,ics=[0,0],ivar=x,end_points=20)
p=line(soln)
p+=plot_vector_field([1/sqrt(1+f(x,y)^2),f(x,y)/sqrt(1+f(x,y)^2)],(x,0,20),(y,0,0.8))
show(p)

more

This didn't work at all until recently, so I'm a little surprised. See Trac 11948 for where this happened...

Uh-oh.

sage: erf(i*55.)
1.00000000000000 + 5.64861461531350e1311*I
sage: import mpmath
sage: mpmath.erf(i*55)
mpc(real='0.0', imag='5.648614615313498e+1311')


It's pretty reliably exactly one off in the real direction for pure imaginary input. That would be a problem. Some other testing reveals that it seems to be fine for other input. This seems, on first glance, to be a bug in the upstream system (Pari) that we still use for evaluating this function.

See Trac 1173 and Trac 13050 for switching this evaluation to mpmath.

more

Also, I'm seeing that the mpmath.erf and erf commands do not give the same imaginary part on some inputs. So, it's off by more than 1.

sage: erf(i*1.42)

1.00000000000000 + 4.03986343036907*I

sage: import mpmath
sage: mpmath.erf(i*1.42)

mpc(real='0.0', imag='3.8217653554366318')

more

Odd. That is really close to the bad spot. I'll add this to the examples.

( 2012-07-02 04:58:35 -0500 )edit

Thanks for all the helpful comments. I now vaguely remember having been unable to plot the error function in Sage before.

Plotting as a vector field works:

maxima.eval('plotdf(-2xy+1,[trajectory_at,0,0],[x,-3,3],[y,-3,3])')

more

You can also use plot_slope_field or plot_vector_field, if you want to stay with Sage commands. See my edit above.

( 2012-07-08 02:15:40 -0500 )edit

Thanks calc314 for the link to Wolfram Alpha. Quite a goldmine, there. Especially liked the 3d (2d real, 1d imaginary) plot of erf. Maybe someday I'll learn how to make a similar plot for the solution of y'+2xy =1.

Your code for the numerical approximation and the slope fields is very interesting. I assume you used 1-xy as a simplification: it doesn't change the shape of the curve in the maxima code plots much. But when I plug in 1-2xy into your code, I have trouble with the y axis coordinates. If I plug 1-2x*y into your code, what is the best way to set the actual plot y coordinates back to say -0.8 to 0.8 and the x coordinates from -10 to 20? Thanks again, Bob

more

I see that I did make that simplification. You can get some improvement by cutting the stepsize in the RK4 algorithm down. This can by done:

y=function('y',x)
f(x,y)=1-2*x*y
soln=desolve_rk4(f(x,y),y,ics=[0,0],ivar=x,end_points=16,stepsize=0.001)
p=line(soln)
p+=plot_vector_field([1/sqrt(1+f(x,y)^2),f(x,y)/sqrt(1+f(x,y)^2)],(x,0,20),(y,0,0.8))
show(p)


You can see that near $x=16$ the approximation starts falling apart. With this sensitivity, I recommend a more robust solver. To use desolve_odeint, you need to have a system. So, below I use $dx/dt=1, dy/dt=1-2xy$ as the system.

var('x,y')
g(x,y)=1-2*x*y
f=[SR(1),g(x,y)]
soln=desolve_odeint(f,[0,0],srange(0,20,0.05),[x,y])
p=line(soln)
p+=plot_vector_field([1/sqrt(1+g(x,y)^2),g(x,y)/sqrt(1+g(x,y)^2)],(x,0,20),(y,0,0.8))
show(p)


(The SR(1) ensures that Sage knows 1 is a symbolic function. desolve_odeint appears to need this.)

To get your plot to start at $x=-10$, you want to start at (0,0) and run back to $x=-10$. Then, combine this with your other plot.

var('x,y')
g(x,y)=1-2*x*y
f=[SR(1),g(x,y)]
soln1=desolve_odeint(f,[0,0],srange(0,-10,-0.05),[x,y])
soln2=desolve_odeint(f,[0,0],srange(0,20,0.05),[x,y])
p=line(soln1)+line(soln2)
p+=plot_vector_field([1/sqrt(1+g(x,y)^2),g(x,y)/sqrt(1+g(x,y)^2)],(x,-10,20),(y,-0.8,0.8))
show(p)


Also, for fun, you might like playing with pplane at http://math.rice.edu/%7Edfield/dfpp.html. This has a very good solver and is handy for some uses. And, it will solve forwards and backwards. (The original pplane program runs in matlab, but this link has a version that runs in java in your browser.)

more

Awesome!

( 2012-07-10 05:40:04 -0500 )edit

Thanks Calc314 for the post up above answering my questions about plotting negative and the link to Rice.edu. I've used pplane with mathlab, but I no longer have mathlab so good to know pplane can be used independently. Bob

more