# sample maple lab 4 Numerical Methods #{Answer Check in Maple}. # Check the exact symbolic solution de:=diff(y(t),t)=0.02225 *y(t) - 0.0003*y(t)^2; ic:= y(0)=25; dsolve({de,ic},y(t)); #{Hand Check for Euler}. # One step. h=2.5 x0 = 0 y0 = 25 f(x,y) = 0.02225 y - 0.0003 y^2 y1 = y0 + h f(x0,y0) = 25 + 2.5 (0.02225 (25) - 0.0003 (25)^2) = 25.921875 Dots[1] = [0, 25], Dots[2] = [2.500000000, 25.92187500]. Answer checks. #{Part III}. Include an appendix of the computer code used. # Now for the Euler code to make the dot table, error percentages and plot. # Euler. Group 1, initialize. f:=(x,y)->0.02225 *y - 0.0003*y^2; x0:=0:y0:=25:Dots:=[x0,y0]:n:=100:h:=250/n: # Group 2, repeat n times. Euler's method for i from 1 to n do Y:=y0+h*f(x0,y0); x0:=x0+h:y0:=Y:Dots:=Dots,[evalf(x0),evalf(y0)]; od: # Group 3, display relevant dots and plot. Exact:=x->2225/(30+59*exp(-89 *x/4000)); P:=unapply(evalf(100*abs(exact-approx)/abs(exact)),(exact,approx)): m:=n/5:X:=[seq(1+m*j,j=0..n/m)]: # List of relevant indices print("Dots"),seq(Dots[k],k=X); print("Exact"),seq(Exact(Dots[k][1]),k=X); print("Error"),seq(P(Exact(Dots[k][1]),Dots[k][2]),k=X); #plot([Dots]); #{Hand Check for Heun}. #One step. h=2.5 x0 = 0 y0 = 25 f(x,y) = 0.02225 y - 0.0003 y^2 y1 = y0 + h f(x0,y0) = 25 + 2.5 (0.02225 (25) - 0.0003 (25)^2) = 25.921875 y2 = y0 + h(f(x0,y0)+f(x0+h,y1))/2 = 25 + 2.5 (0.02225 (25) - 0.0003 (25)^2)/2 + 2.5 (0.02225 (25.921875) - 0.0003 (25.921875)^2)/2 = 25.92991080 Dots[1] = [0, 25], Dots[2] = [2.500000000, 25.92991080]. Answer checks. #{Part III}. Include an appendix of the computer code used. # Now for the Heun code to make the dot table, error percentages and plot. # Heun. Group 1, initialize. f:=(x,y)->0.02225 *y - 0.0003*y^2; x0:=0:y0:=25:Dots:=[x0,y0]:n:=100:h:=250/n: # Group 2, repeat n times. Heun's method for i from 1 to n do Y1:=y0+h*f(x0,y0); Y:=y0+h*(f(x0,y0)+f(x0+h,Y1))/2; x0:=x0+h:y0:=Y:Dots:=Dots,[evalf(x0),evalf(y0)]; od: # Group 3, display relevant dots and plot. Exact:=x->2225/(30+59*exp(-89 *x/4000)); P:=unapply(evalf(100*abs(exact-approx)/abs(exact)),(exact,approx)): m:=n/5:X:=[seq(1+m*j,j=0..n/m)]: # List of relevant indices print("Dots"),seq(Dots[k],k=X); print("Exact"),seq(Exact(Dots[k][1]),k=X); print("Error"),seq(P(Exact(Dots[k][1]),Dots[k][2]),k=X); #plot([Dots]); #{Part III}. Include an appendix of the computer code used. # Now for the RK4 code to make the dot table, error percentages and plot. # RK4. Group 1, initialize. f:=(x,y)->0.02225 *y - 0.0003*y^2; x0:=0:y0:=25:Dots:=[x0,y0]:n:=100:h:=250/n: # Group 2, repeat n times. RK4 method. for i from 1 to n do k1:=h*f(x0,y0): k2:=h*f(x0+h/2,y0+k1/2): k3:=h*f(x0+h/2,y0+k2/2): k4:=h*f(x0+h,y0+k3): Y:=y0+(k1+2*k2+2*k3+k4)/6: x0:=x0+h:y0:=Y:Dots:=Dots,[evalf(x0),evalf(y0)]; od: # Group 3, display relevant dots and plot. Exact:=x->2225/(30+59*exp(-89 *x/4000)); P:=unapply(evalf(100*abs(exact-approx)/abs(exact)),(exact,approx)): m:=n/5:X:=[seq(1+m*j,j=0..n/m)]: # List of relevant indices print("Dots"),seq(Dots[k],k=X); print("Exact"),seq(Exact(Dots[k][1]),k=X); print("Error"),seq(P(Exact(Dots[k][1]),Dots[k][2]),k=X); plot([Dots]); # Warning: These snips of code made for y'=1-x-y, y(0)=3. # Code computes approx values for y(0.1) to y(1). # 'Dots' is the list of dots for connect-the-dots graphics. # ======================================== # Euler. Group 1, initialize. f:=(x,y)->1-x-y: x0:=0:y0:=3:h:=0.1:Dots:=[x0,y0]:n:=200: # Group 2, repeat n times. Euler's method for i from 1 to n do Y:=y0+h*f(x0,y0); x0:=x0+h:y0:=Y:Dots:=Dots,[x0,y0]; od: # Group 3, display relevant dots and plot. Exact:=x->2-x+exp(-x); P:=unapply(evalf(100*abs(exact-approx)/abs(exact)),(exact,approx)): m:=40:X:=[seq(1+m*j,j=0..n/m)]: # List of relevant indices print("Dots"),seq(Dots[k],k=X); print("Exact"),seq(Exact(Dots[k][1]),k=X); print("Error"),seq(P(Exact(Dots[k][1]),Dots[k][2]),k=X); plot([Dots]); # ======================================== # Heun. Group 1, initialize. f:=(x,y)->1-x-y: x0:=0:y0:=3:h:=0.1:Dots:=[x0,y0]:n:=200: # Group 2, repeat n times. Heun method. for i from 1 to n do Y1:=y0+h*f(x0,y0); Y:=y0+h*(f(x0,y0)+f(x0+h,Y1))/2: x0:=x0+h:y0:=Y:Dots:=Dots,[x0,y0]; od: # Group 3, display relevant dots and plot. Dots[1],Dots[2],seq(Dots[1+40*j],j=1..n/40); plot([Dots]); # ======================================== # RK4. Group 1, initialize. f:=(x,y)->1-x-y: x0:=0:y0:=3:h:=0.1:Dots:=[x0,y0]:n:=100: # Group 2, repeat n times. RK4 method. for i from 1 to n do k1:=h*f(x0,y0): k2:=h*f(x0+h/2,y0+k1/2): k3:=h*f(x0+h/2,y0+k2/2): k4:=h*f(x0+h,y0+k3): Y:=y0+(k1+2*k2+2*k3+k4)/6: x0:=x0+h:y0:=Y:Dots:=Dots,[x0,y0]; od: # Group 3, display some dots and plot. Dots[1],Dots[2],Dots[101]; plot([Dots]);