Plot solution for y' + 2xy = 1

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)

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

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?

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

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.