ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 11 Jul 2012 09:24:39 -0500Plot solution for y' + 2xy = 1http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/f = desolve(diff(y,x) + 2*x*y - 1, y, ics=[0,0]); f
plot(f) # error message ... unable to simplify to float approximation.
Tried plotting real_part of solution: -1/2*I*sqrt(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.Sun, 01 Jul 2012 15:29:48 -0500http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/Answer by calc314 for <p>f = desolve(diff(y,x) + 2<em>x</em>y - 1, y, ics=[0,0]); f
plot(f) # error message ... unable to simplify to float approximation.</p>
<p>Tried plotting real_part of solution: -1/2<em>I</em>sqrt(pi)<em>e^(-x^2)</em>erf(I*x),
but get same error message.</p>
<p>Tried using list_plot, but error about symbolic expression.
Haven't been able to get implicit_plot or sol.simplify_full to work.</p>
<p>Thanks for any hints.</p>
http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13781#post-id-13781Your 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](http://m.wolframalpha.com/input/?i=plot+-1%2F2*I*sqrt%28pi%29*e%5E%28-x%5E2%29*erf%28I*x%29+from+0+to+10&x=0&y=0).)
Anyone have ideas on this?Mon, 02 Jul 2012 01:29:31 -0500http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13781#post-id-13781Comment by kcrisman for <p>Your solution is <code>-1/2*I*sqrt(pi)*e^(-x^2)*erf(I*x)</code> which Sage does produce correctly. But, Sage appears to be having trouble with <code>erf(I*x)</code>. This should be a pure imaginary result, but Sage is not giving that for $x>1.42$.</p>
<p>Wolfram alpha plots your function nicely and also handles <code>erf(I*x)</code> well. (See <a href="http://m.wolframalpha.com/input/?i=plot+-1%2F2*I*sqrt%28pi%29*e%5E%28-x%5E2%29*erf%28I*x%29+from+0+to+10&x=0&y=0">this link</a>.)</p>
<p>Anyone have ideas on this?</p>
http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?comment=19459#post-id-19459Thank you very much for finding this!Mon, 02 Jul 2012 03:59:24 -0500http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?comment=19459#post-id-19459Comment by calc314 for <p>Your solution is <code>-1/2*I*sqrt(pi)*e^(-x^2)*erf(I*x)</code> which Sage does produce correctly. But, Sage appears to be having trouble with <code>erf(I*x)</code>. This should be a pure imaginary result, but Sage is not giving that for $x>1.42$.</p>
<p>Wolfram alpha plots your function nicely and also handles <code>erf(I*x)</code> well. (See <a href="http://m.wolframalpha.com/input/?i=plot+-1%2F2*I*sqrt%28pi%29*e%5E%28-x%5E2%29*erf%28I*x%29+from+0+to+10&x=0&y=0">this link</a>.)</p>
<p>Anyone have ideas on this?</p>
http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?comment=19458#post-id-19458Glad to help!Mon, 02 Jul 2012 04:42:08 -0500http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?comment=19458#post-id-19458Answer by bbtp for <p>f = desolve(diff(y,x) + 2<em>x</em>y - 1, y, ics=[0,0]); f
plot(f) # error message ... unable to simplify to float approximation.</p>
<p>Tried plotting real_part of solution: -1/2<em>I</em>sqrt(pi)<em>e^(-x^2)</em>erf(I*x),
but get same error message.</p>
<p>Tried using list_plot, but error about symbolic expression.
Haven't been able to get implicit_plot or sol.simplify_full to work.</p>
<p>Thanks for any hints.</p>
http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13810#post-id-13810Thanks 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-x*y as a simplification: it doesn't change the shape of the curve in the maxima code plots much.
But when I plug in 1-2*x*y into your code, I have trouble with the y axis coordinates.
If I plug 1-2*x*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,
BobMon, 09 Jul 2012 15:17:49 -0500http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13810#post-id-13810Answer by kcrisman for <p>f = desolve(diff(y,x) + 2<em>x</em>y - 1, y, ics=[0,0]); f
plot(f) # error message ... unable to simplify to float approximation.</p>
<p>Tried plotting real_part of solution: -1/2<em>I</em>sqrt(pi)<em>e^(-x^2)</em>erf(I*x),
but get same error message.</p>
<p>Tried using list_plot, but error about symbolic expression.
Haven't been able to get implicit_plot or sol.simplify_full to work.</p>
<p>Thanks for any hints.</p>
http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13783#post-id-13783This didn't work at all until recently, so I'm a little surprised. See [Trac 11948](http://trac.sagemath.org/sage_trac/ticket/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](http://trac.sagemath.org/sage_trac/ticket/1173) and [Trac 13050](http://trac.sagemath.org/sage_trac/ticket/13050) for switching this evaluation to mpmath.Mon, 02 Jul 2012 03:59:02 -0500http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13783#post-id-13783Answer by bbtp for <p>f = desolve(diff(y,x) + 2<em>x</em>y - 1, y, ics=[0,0]); f
plot(f) # error message ... unable to simplify to float approximation.</p>
<p>Tried plotting real_part of solution: -1/2<em>I</em>sqrt(pi)<em>e^(-x^2)</em>erf(I*x),
but get same error message.</p>
<p>Tried using list_plot, but error about symbolic expression.
Haven't been able to get implicit_plot or sol.simplify_full to work.</p>
<p>Thanks for any hints.</p>
http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13813#post-id-13813Thanks 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.
BobWed, 11 Jul 2012 09:24:39 -0500http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13813#post-id-13813Answer by calc314 for <p>f = desolve(diff(y,x) + 2<em>x</em>y - 1, y, ics=[0,0]); f
plot(f) # error message ... unable to simplify to float approximation.</p>
<p>Tried plotting real_part of solution: -1/2<em>I</em>sqrt(pi)<em>e^(-x^2)</em>erf(I*x),
but get same error message.</p>
<p>Tried using list_plot, but error about symbolic expression.
Haven't been able to get implicit_plot or sol.simplify_full to work.</p>
<p>Thanks for any hints.</p>
http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13785#post-id-13785Also, 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')Mon, 02 Jul 2012 04:41:57 -0500http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13785#post-id-13785Comment by kcrisman for <p>Also, I'm seeing that the <code>mpmath.erf</code> and <code>erf</code> commands do not give the same imaginary part on some inputs. So, it's off by more than 1.</p>
<pre><code>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')
</code></pre>
http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?comment=19457#post-id-19457Odd. That is really close to the bad spot. I'll add this to the examples.Mon, 02 Jul 2012 04:58:35 -0500http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?comment=19457#post-id-19457Answer by bbtp for <p>f = desolve(diff(y,x) + 2<em>x</em>y - 1, y, ics=[0,0]); f
plot(f) # error message ... unable to simplify to float approximation.</p>
<p>Tried plotting real_part of solution: -1/2<em>I</em>sqrt(pi)<em>e^(-x^2)</em>erf(I*x),
but get same error message.</p>
<p>Tried using list_plot, but error about symbolic expression.
Haven't been able to get implicit_plot or sol.simplify_full to work.</p>
<p>Thanks for any hints.</p>
http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13793#post-id-13793
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(-2*x*y+1,[trajectory_at,0,0],[x,-3,3],[y,-3,3])')
Sat, 07 Jul 2012 17:14:51 -0500http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13793#post-id-13793Comment by calc314 for <p>Thanks for all the helpful comments.
I now vaguely remember having been unable to plot the error function in Sage before.</p>
<p>Plotting as a vector field works:</p>
<p>maxima.eval('load("plotdf")')</p>
<p>maxima.eval('plotdf(-2<em>x</em>y+1,[trajectory_at,0,0],[x,-3,3],[y,-3,3])')</p>
http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?comment=19427#post-id-19427You can also use `plot_slope_field` or `plot_vector_field`, if you want to stay with Sage commands. See my edit above.Sun, 08 Jul 2012 02:15:40 -0500http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?comment=19427#post-id-19427Answer by calc314 for <p>f = desolve(diff(y,x) + 2<em>x</em>y - 1, y, ics=[0,0]); f
plot(f) # error message ... unable to simplify to float approximation.</p>
<p>Tried plotting real_part of solution: -1/2<em>I</em>sqrt(pi)<em>e^(-x^2)</em>erf(I*x),
but get same error message.</p>
<p>Tried using list_plot, but error about symbolic expression.
Haven't been able to get implicit_plot or sol.simplify_full to work.</p>
<p>Thanks for any hints.</p>
http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13811#post-id-13811I 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](/upfiles/13418862593629736.png)
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](/upfiles/13418864966062555.png)
(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](/upfiles/13418896778217787.png)
Also, for fun, you might like playing with `pplane` at [http://math.rice.edu/%7Edfield/dfpp.html](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.)
Mon, 09 Jul 2012 16:19:15 -0500http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13811#post-id-13811Comment by kcrisman for <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:</p>
<pre><code>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)
</code></pre>
<p><img alt="image description" src="/upfiles/13418862593629736.png"/></p>
<p>You can see that near $x=16$ the approximation starts falling apart. With this sensitivity, I recommend a more robust solver. To use <code>desolve_odeint</code>, you need to have a system. So, below I use $dx/dt=1, dy/dt=1-2xy$ as the system.</p>
<pre><code>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)
</code></pre>
<p><img alt="image description" src="/upfiles/13418864966062555.png"/></p>
<p>(The <code>SR(1)</code> ensures that Sage knows <code>1</code> is a symbolic function. <code>desolve_odeint</code> appears to need this.)</p>
<p>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.</p>
<pre><code>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)
</code></pre>
<p><img alt="image description" src="/upfiles/13418896778217787.png"/></p>
<p>Also, for fun, you might like playing with <code>pplane</code> at <a href="http://math.rice.edu/%7Edfield/dfpp.html">http://math.rice.edu/%7Edfield/dfpp.html</a>. 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.)</p>
http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?comment=19426#post-id-19426Awesome! Tue, 10 Jul 2012 05:40:04 -0500http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?comment=19426#post-id-19426Answer by calc314 for <p>f = desolve(diff(y,x) + 2<em>x</em>y - 1, y, ics=[0,0]); f
plot(f) # error message ... unable to simplify to float approximation.</p>
<p>Tried plotting real_part of solution: -1/2<em>I</em>sqrt(pi)<em>e^(-x^2)</em>erf(I*x),
but get same error message.</p>
<p>Tried using list_plot, but error about symbolic expression.
Haven't been able to get implicit_plot or sol.simplify_full to work.</p>
<p>Thanks for any hints.</p>
http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13782#post-id-13782You 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)
Mon, 02 Jul 2012 01:44:27 -0500http://ask.sagemath.org/question/9121/plot-solution-for-y-2xy-1/?answer=13782#post-id-13782