Ask Your Question
0

Plot solution for y' + 2xy = 1

asked 2012-07-01 15:29:48 -0500

bbtp gravatar image

updated 2015-01-24 06:55:06 -0500

FrédéricC gravatar image

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 flag offensive close merge delete

8 answers

Sort by » oldest newest most voted
3

answered 2012-07-02 01:44:27 -0500

calc314 gravatar image

updated 2012-07-08 02:15:00 -0500

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)
edit flag offensive delete link more
2

answered 2012-07-09 16:19:15 -0500

calc314 gravatar image

updated 2012-07-09 17:08:22 -0500

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)

image description

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)

image description

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

image description

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

edit flag offensive delete link more

Comments

Awesome!

kcrisman gravatar imagekcrisman ( 2012-07-10 05:40:04 -0500 )edit
1

answered 2012-07-02 01:29:31 -0500

calc314 gravatar image

updated 2012-07-02 01:32:40 -0500

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?

edit flag offensive delete link more

Comments

Thank you very much for finding this!

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

Glad to help!

calc314 gravatar imagecalc314 ( 2012-07-02 04:42:08 -0500 )edit
1

answered 2012-07-02 03:59:02 -0500

kcrisman gravatar image

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.

edit flag offensive delete link more
1

answered 2012-07-02 04:41:57 -0500

calc314 gravatar image

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')
edit flag offensive delete link more

Comments

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

kcrisman gravatar imagekcrisman ( 2012-07-02 04:58:35 -0500 )edit
0

answered 2012-07-09 15:17:49 -0500

bbtp gravatar image

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

edit flag offensive delete link more
0

answered 2012-07-07 17:14:51 -0500

bbtp gravatar image

updated 2012-07-07 17:18:12 -0500

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('load("plotdf")')

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

edit flag offensive delete link more

Comments

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

calc314 gravatar imagecalc314 ( 2012-07-08 02:15:40 -0500 )edit
0

answered 2012-07-11 09:24:39 -0500

bbtp gravatar image

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

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2012-07-01 15:29:48 -0500

Seen: 1,549 times

Last updated: Jul 11 '12