1426.sagews

In order to be able to approximate solution to autonomous DE using Eulers method, we need to specifically tell sage about the dependence between the variables x and y.  If we don't do that, sage will think that expressions like y*(1-y) are just functions of variable y, which, as you have seen in class, completely confuses the eulers_method function.

One way to do that is using the PolynomialRing function.  The following will tell sage that x and y are both generators of the same two variable polynomial ring over Q.  That seems to take care of the eulers_method call for autonomous equations, as long as the right hand side is a polynomial.  Autonomous equations with right hand sides like sin(y) will still fail.  I was not able to find a general solution so far.

In other words, don't worry about the following, the purpose is to allow us to use eulers_method for at least some autonomous equations.

If your equations are not autonomous, it seems like it is enough to just use "y = var('y')", declaring y as a symbolic variable.

x,y = PolynomialRing(QQ,2,"xy").gens()

Start by plotting slope field for

y' = y - x

plot_slope_field(y-x,(x,-4,4),(y,-4,4)).show(aspect_ratio=1)

Now we will save the plot into a variable to use it later:

sf = plot_slope_field(y-x,(x,-4,4),(y,-4,4))

Now we will use Euler's method to draw approximate solution curve.  Again, we will save the plots to variables, so we can later combine them with the slope field we created above.

p1 = list_plot(filter(lambda p: p[1]<4 and p[1]>-4,eulers_method(y-x,-3,-1.2,.01,4,method="none")))

(The expression "lambda p: p[1]<4 and p[1]>-4" is an anonumous function which takes a single argument "p", and returns True if the second component of p is between -4 and 4.  Note that the second component is p[1], the first component would be p[0]. 

The function "filter" takes two variables, a function and a list.  It goes through the list, feeding each element into the function, and keeping the element if the function returns "true", and throwing it away if the function returns "false".

So the constuct above takes the list of points returned by eulers_method, and remove all of the points whose y-coordinate is not between -4 an 4.)

p2 = list_plot(filter(lambda x: x[1]<4 and x[1]>-4,eulers_method(y-x,-3,-2,.01,4,method="none")))

p3 = list_plot(filter(lambda x: x[1]<4 and x[1]>-4,eulers_method(y-x,-3,-2.2,.01,4,method="none")))

Now we will show all the plots we created and saved in variables:

show(sf+p1+p2+p3,aspect_ratio=1)

Now we will write a sage command that will do all the work we did above for us. The procedure "deplot" will have 5 mandatory arguments:

• the expression "f": the right hand side of the normal for of our DE
• the xmin: the lower bound of the viewing window x coordinate
• the xmax: the upper bound of the viewing window x coordinate
• the ymin: the lower bound of the viewing window y coordinate
• the ymax: the upper bound of the viewing window y coordinate

and an optional argument ivlist, containing list of the initial points.

def plotde (f, xmin, xmax, ymin, ymax, ivlist=[]):
    theplot = plot_slope_field(f,(x,xmin,xmax),(y,ymin,ymax))
    for pt in ivlist:
        theplot = theplot + list_plot(filter(lambda x: x[1]<ymax and x[1]>ymin,eulers_method(f,pt[0],pt[1],.01,xmax,method="none")))
    return(theplot)

Test it by plottin the slope field of y'=x+y for x between -3 and 3, y between -3 and 3, together with approximate solution curves with initial conditions

• y(-2)=-2
• y(-2)=1
• y(-3)=1.9

plotde(x+y,-3,3,-3,3,ivlist=[(-2,2),(-2,-1),(-3,1.9)]).show(aspect_ratio=1)

This time y'=4y/x, with initial conditions y(-3)=2, y(-2)=-1 and y(-.1)=2

plotde(4*y/x,-3,3,-3,3,ivlist=[(-3,2),(-2,-1),(-.1,2)]).show(aspect_ratio=1)

Another one, with some rounding error mess (can you explain what's going on?):

plotde(-x/y,-3,3,-3,3,ivlist=[(1,1),(-2,-1)]).show(aspect_ratio=1)

And finally an autonomous equation, the one I couldn't get working in class:

plotde(y*(1-y) ,-3,3,-3,3,ivlist=[(0,.5),(0,1),(0,1.2)]).show(aspect_ratio=1)