4309.sagews

t=var("t");
f(t, y) = (9/10)*y*t - (2/12)*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.

22/27*e^(9/20*t^2) + 5/27

# 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, t0, y0):
    antideriv = integral(f(t,y), t)
    return (y0 + antideriv - antideriv(t=t0)).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 < 11 :  #Computes up to the iterate p_5.
    # how close is each iterate to the actual value at t=2?
    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, 22/27*e^(36/5) + 5/27, difference == 1090.5732154324878) (1, 103/15, 22/27*e^(36/5) + 5/27, difference == 1084.706548765821) (2, 2099/75, 22/27*e^(36/5) + 5/27, difference == 1063.5865487658211) (3, 29503/375, 22/27*e^(36/5) + 5/27, difference == 1012.8985487658211) (4, 318587/1875, 22/27*e^(36/5) + 5/27, difference == 921.66014876582108) (5, 14123267/46875, 22/27*e^(36/5) + 5/27, difference == 790.2768527658211) (6, 107567887/234375, 22/27*e^(36/5) + 5/27, difference == 632.61689756582109) (7, 5095131917/8203125, 22/27*e^(36/5) + 5/27, difference == 470.45237221724972) (8, 31461811009/41015625, 22/27*e^(36/5) + 5/27, difference == 324.50429940353547) (9, 181253660741/205078125, 22/27*e^(36/5) + 5/27, difference == 207.74584115256403) (10, 4962344421053/5126953125, 22/27*e^(36/5) + 5/27, difference == 123.67975121186464) (11, 26222277058993/25634765625, 22/27*e^(36/5) + 5/27, difference == 68.654674159770479)

G=plot((p[0], p[1], p[2], p[3], p[4], p[5],p[6],p[7],p[8],p[9],p[10],p[11]), 0, 4, rgbcolor=hue(1));
G+=P_g;
G.show(figsize=5)