MA 531 HW 7 Sample Code.sagews

#Function I wrote for 7.10.x.  Adapt it to work on 7.10.8.

#Code for whole table, more accurate than book's.
#Note: yP need to be isolated to y'= expression.
#All derivatives also need to be taken and put in terms of y'', y''', etc. and set as y2P,y3P,y4P, and y5P.
def PrintTaylor5(a,b,n,y0,yP,y2P,y3P,y4P,y5P):
    h=(b-a)/n
    xi=a
    yi=y0
    print "y("+str(a.n())+")=",y0.n()
    for i in [1,2,..,n]:
        newyi=yi+h*yP(xi,yi)+h^2/2*y2P(xi,yi)+h^3/(factorial(3))*y3P(xi,yi)+h^4/(factorial(4))*y4P(xi,yi)+h^5/(factorial(5))*y5P(xi,yi)
        xi=xi+h
        yi=newyi
        print "y("+str(xi.n())+")=",yi.n()
a=0
b=1
n=5
y0=7
var('x y')
yP(x,y)=2*y+3*sin(x)+e^(x)
y2P(x,y)=2*yP+3*cos(x)+e^(x)
y3P(x,y)=2*y2P-3*sin(x)+e^(x)
y4P(x,y)=2*y3P-3*cos(x)+e^(x)
y5P(x,y)=2*y4P+3*sin(x)+e^(x)
PrintTaylor5(a,b,n,y0,yP,y2P,y3P,y4P,y5P)

#7.10.14. [Edited.] Use the Sage library method to solve y'''(x)=cos(2y)+sin(y')-e^(y'')+x^5, with initial conditions y(1)=4, y'(1)=7, and y''(1)=-2 to find the value of y(1.5) with n=50 steps.

from sage.calculus.desolvers import desolve_system_rk4

#Note there should be a lower case sage in the import command (typo in book)