4342.sagews

t=var("t");
f(t,y) = .6*y*t-.8*t ;
y = function("y", t) ; #defines a variable t and a function y which depends on t

g=desolve(diff(y,t)-f(t,y) , y, ics=[0,1]); g  

#g is the solution to the ODE dy/dt - y - t = 0 with initial condition t=0, y=1. 

print(g); #prints the solution

P_g=plot(g, 0, 4);P_g.show(figsize=4) #plots the solution between t=0 and t=4.

-1/3*e^(3/10*t^2) + 4/3

# Compute the next Picard iteration.
# f(t,y) should be the RHS of the DE
# y(t) should be the function to Picard-iterate
# (t0,y0) is the initial value

def PicardIteration (f, y, 0, 1):
    antideriv = integral(f(t,y), t)
    return (1 + antideriv - antideriv(t=0)).simplify()

# define our first iteration to be constantly equal to the initial value
p0(t) = 1           # initial Picard iteration, identically = 1
p = [p0]            # list of Picard iterations (p for Picard)


while len(p)-2 < 6 :  #Computes up to the iterate p_5.
    # how close is each iterate to the actual value at t=4?
    print(len(p)-1, (p[-1])(4), g(4), var('difference') == float( abs((p[-1])(4) - g(4))))
    p.append(PicardIteration(f, p[-1], 0, 1))

(0, 1, -1/3*e^(24/5) + 4/3, difference == 40.170139172911611) (1, -0.6, -1/3*e^(24/5) + 4/3, difference == 38.57013917291161) (2, -4.44, -1/3*e^(24/5) + 4/3, difference == 34.730139172911613) (3, -10.584, -1/3*e^(24/5) + 4/3, difference == 28.586139172911615) (4, -17.9568, -1/3*e^(24/5) + 4/3, difference == 21.213339172911613) (5, -25.034688, -1/3*e^(24/5) + 4/3, difference == 14.135451172911615) (6, -30.6969984, -1/3*e^(24/5) + 4/3, difference == 8.4731407729116128)

G=plot((p[0], p[1], p[2], p[3], p[4], p[5]), 0, 2, rgbcolor=hue(1));
G+=P_g;
G.show(figsize=5)