> # This worksheet demonstrates the use of the "Share" library routine # firsteuler. It graphs the Euler method solution for h = .1,.05,.025, # and .01 along with the actual solution from dsolve. It then plots all # of these together to show a comparision. -------------------------------------------------------------------------------- # # Art Belmonte # Mon, 29/Jan/96 # Math 308-509 # Chapter 8: Numerical Methods (for solving ODEs) # Euler's Method [ T-389/Eq (6) ] -------------------------------------------------------------------------------- > with(share):\ readshare(ODE,plots): See ?share and ?share,contents for information about the share library > ?firsteuler -------------------------------------------------------------------------------- # T-390/Table 8.1.1: Column 1, via raw dump; graph of solution (all # plots omitted to save space). # LEGEND: # a - left endpoint of interval for independent variable # b - right endpoint of interval for independent variable # f - derivative function in D.E.; from dy/dt = f(t, y) # h - step size # i - list of initial values: [ t0, y0 ] = [ indep_var, dep_var ] # n - number of steps: (b-a)/h, as a RATIONAL! > unassign('t', 'y');\ a:=0; b:=1; h:=0.1; i:=[0, 1];\ n:=convert((b-a)/h, rational);\ f:=(t, y)->1-t+4*y;\ deq:={diff(y(t), t) = f(t, y(t))};\ eulerpts:= firsteuler(f, i, h, n);\ list1:=convert(eulerpts,listlist):\ plot(list1); a := 0 b := 1 h := .1 i := [0, 1] n := 10 f := (t,y) -> 1 - t + 4 y d deq := {---- y(t) = 1 - t + 4 y(t)} dt eulerpts := array(0 .. 10,, [ 0 = [0, 1.] 1 = [.1, 1.5] 2 = [.2, 2.19] 3 = [.3, 3.146] 4 = [.4, 4.4744] 5 = [.5, 6.32416] 6 = [.6, 8.903824] 7 = [.7, 12.5053536] 8 = [.8, 17.53749504] 9 = [.9, 24.57249306] 10 = [1.0, 34.41149028] ]) -------------------------------------------------------------------------------- # T-390/Table 8.1.1: Column 2, via judicious paring (every 2nd # item); graph of solution. > unassign('t', 'y');\ a:=0; b:=1; h:=0.05; i:=[0, 1];\ n:=convert((b-a)/h, rational);\ f:=(t, y)->1-t+4*y;\ deq:={diff(y(t), t) = f(t, y(t))};\ eulerpts:=firsteuler(f, i, h, n):\ for i from 0 to n by 2 do\ eulerpts[i]\ od;\ list2 := convert(eulerpts,listlist):\ plot(list2); a := 0 b := 1 h := .05 i := [0, 1] n := 20 f := (t,y) -> 1 - t + 4 y d deq := {---- y(t) = 1 - t + 4 y(t)} dt [0, 1.] [.10, 1.5475] [.20, 2.32490000] [.30, 3.433356000] [.40, 5.018532640] [.50, 7.290187002] [.60, 10.55036928] [.70, 15.23403177] [.80, 21.96750574] [.90, 31.65270827] [1.00, 45.58839992] -------------------------------------------------------------------------------- # T-390/Table 8.1.1: Column 3, via judicious paring (every 4th item); # graph of solution. > unassign('t', 'y');\ a:=0; b:=1; h:=0.025; i:=[0, 1];\ n:=convert((b-a)/h, rational);\ f:=(t, y)->1-t+4*y;\ deq:={diff(y(t), t) = f(t, y(t))};\ eulerpts:=firsteuler(f, i, h, n):\ for i from 0 to n by 4 do\ eulerpts[i]\ od;\ list3:=convert(eulerpts,listlist):\ plot(list3); a := 0 b := 1 h := .025 i := [0, 1] n := 40 f := (t,y) -> 1 - t + 4 y d deq := {---- y(t) = 1 - t + 4 y(t)} dt [0, 1.] [.100, 1.576118750] [.200, 2.408011713] [.300, 3.614383699] [.400, 5.369030425] [.500, 7.926406197] [.600, 11.65905757] [.700, 17.11242994] [.800, 25.08510991] [.900, 36.74630817] [1.000, 53.80786605] -------------------------------------------------------------------------------- # T-390/Table 8.1.1: Column 4, via judicious paring (every 10th # item); graph of solution. > unassign('t', 'y');\ a:=0; b:=1; h:=0.01; i:=[0, 1];\ n:=convert((b-a)/h, rational);\ f:=(t, y)->1-t+4*y;\ deq:={diff(y(t), t) = f(t, y(t))};\ eulerpts:=firsteuler(f, i, h, n):\ for i from 0 to n by 10 do\ eulerpts[i]\ od;\ list4:=convert(eulerpts,listlist):\ plot(list4); a := 0 b := 1 h := .01 i := [0, 1] n := 100 f := (t,y) -> 1 - t + 4 y d deq := {---- y(t) = 1 - t + 4 y(t)} dt [0, 1.] [.10, 1.595290089] [.20, 2.464458735] [.30, 3.739034547] [.40, 5.613712001] [.50, 8.376686480] [.60, 12.45455756] [.70, 18.47879679] [.80, 27.38413639] [.90, 40.55420836] [1.00, 60.03712599] -------------------------------------------------------------------------------- # T-390/Table 8.1.1: Column 5, exact solution (via dsolve); graph of # solution. > unassign('t', 'y');\ f:=(t, y)->1-t+4*y;\ deq:={diff(y(t), t) = f(t, y(t))};\ IC:={y(0)=1}; IVP:=deq union IC;\ sol:=dsolve(IVP, y(t));\ y:=unapply(subs(sol, y(t)), t);\ check:=IVP;\ for i from 0.0 to 1.0 by 0.1 do\ [i, y(i)]\ od;\ plot(y(t), t=0..1); f := (t,y) -> 1 - t + 4 y d deq := {---- y(t) = 1 - t + 4 y(t)} dt IC := {y(0) = 1} d IVP := {y(0) = 1, ---- y(t) = 1 - t + 4 y(t)} dt 19 sol := y(t) = - 3/16 + 1/4 t + ---- exp(4 t) 16 y := t -> - 3/16 + 1/4 t + 19/16 exp(4 t) check := {1 = 1, 1/4 + 19/4 exp(4 t) = 1/4 + 19/4 exp(4 t)} [0, 1] [.1, 1.609041829] [.2, 2.505329852] [.3, 3.830138846] [.4, 5.794226004] [.5, 8.712004118] [.6, 13.05252195] [.7, 19.51551804] [.8, 29.14487961] [.9, 43.49790340] [1.0, 64.89780316] -------------------------------------------------------------------------------- # The plots from the top down are # 1) actual solution # 2) euler with h=0.01 # 3) euler with h=0.025 # 4) euler with h=0.05 # 5) euler with h=0.1 > plot( {list1,list2,list3,list4,y(t)},t=0..1); >