637.sagews

#adams-bashforth two-step method

var('t','y')
f(t,y) = -3*y
start = 0.0
end = 1.0
solution(t) = float(e)^(-3*t)

numOfIterations = 20

h = (end-start)/numOfIterations

y00 = 1.0 #initial value

t00 = start;

soln = [[t00,y00]]

#trapezoid rule to get started
y01 = y00 + h/2*(f(t00,y00) + f(t00+h,y00 + h*f(t00,y00)))
t01 = t00 + h
soln.append([t01,y01])

for i in range (1,numOfIterations):
    y02 = y01 + h*(3/2*f(t01,y01) - 1/2*f(t00,y00))
    t02 = t01 + h
    t01, t00 = t02, t01
    y01, y00 = y02, y01
    soln.append([t02,y02])
    print t02, y02, solution(t02)

p = list_plot(soln, rgbcolor=(1,0,0))
p1 = plot(solution(x),0,1)
(p + p1).show()

0.100000000000000 0.742468750000000 0.740818220681718 0.150000000000000 0.640007031250000 0.637628151621773 0.200000000000000 0.551690605468750 0.548811636094026 0.250000000000000 0.475560746582031 0.472366552741015 0.300000000000000 0.409936374011231 0.406569659740599 0.350000000000000 0.353367745852356 0.349937749111155 0.400000000000000 0.304605231086418 0.301194211912202 0.450000000000000 0.262571635030901 0.259240260645892 0.500000000000000 0.226338409480429 0.223130160148430 0.550000000000000 0.195105139974650 0.192049908620754 0.600000000000000 0.168181864191386 0.165298888221587 0.650000000000000 0.144973830246423 0.142274071586514 0.700000000000000 0.124968358255332 0.122456428252982 0.750000000000000 0.107723514916364 0.105399224561864 0.800000000000000 0.0928583509293320 0.0907179532894125 0.850000000000000 0.0800444855889596 0.0780816660011531 0.900000000000000 0.0689988526511436 0.0672055127397497 0.950000000000000 0.0594774472238082 0.0578443208748384 1.00000000000000 0.0512699355472871 0.0497870683678639

#unstable Adams-Bashforth (like)

var('t','y')
f(t,y) = -3*y
start = 0.0
end = 1.0
solution(t) = float(e)^(-3*t)

numOfIterations = 20

h = (end-start)/numOfIterations

y00 = 1.0 #initial value

t00 = start;

soln = [[t00,y00]]

#trapezoid rule to get started
y01 = y00 + h/2*(f(t00,y00) + f(t00+h,y00 + h*f(t00,y00)))
t01 = t00 + h
soln.append([t01,y01])

for i in range (1,numOfIterations):
    y02 = -y01 + 2*y00 + h*(5/2*f(t01,y01) + 1/2*f(t00,y00))
    t02 = t01 + h
    t01, t00 = t02, t01
    y01, y00 = y02, y01
    soln.append([t02,y02])
    print t02, y02, solution(t02)

p = list_plot(soln, rgbcolor=(1,0,0), plotjoined=True)
p2 = list_plot(soln, rgbcolor=(1,0,0))
p1 = plot(solution(x),0,1)
(p + p1+p2).show(ymax=5,ymin=-5)

0.100000000000000 0.740781250000000 0.740818220681718 0.150000000000000 0.639332031250000 0.637628151621773 0.200000000000000 0.546922363281250 0.548811636094026 0.250000000000000 0.478695910644531 0.472366552741015 0.300000000000000 0.394618672180176 0.406569659740599 0.350000000000000 0.378888953742981 0.349937749111155 0.400000000000000 0.238668632550239 0.301194211912202 0.450000000000000 0.401191866198661 0.259240260645892 0.500000000000000 -0.0922016983639490 0.223130160148430 0.550000000000000 0.899071677682852 0.192049908620754 0.600000000000000 -1.41371182616452 0.165298888221587 0.650000000000000 3.67456674051571 0.142274071586514 0.700000000000000 -7.77392453357581 0.122456428252982 0.750000000000000 17.7626872091595 0.105399224561864 0.800000000000000 -39.3884996397277 0.0907179532894125 0.850000000000000 88.3523598822576 0.0780816660011531 0.900000000000000 -197.307356644580 0.0672055127397497 0.950000000000000 441.375908159643 0.0578443208748384 1.00000000000000 -986.708535260326 0.0497870683678639