#Kafung Mok #HW 16 #March 30, 2014. ########################################### Class Notes March 27rd ################################################ #March 24, 2014, multi-step methods for solving ODEs, implicit Runge-Kutta Help:=proc(): print(`AAB(f,x,y,x0,y0,h,x1,L), FAB(d) , GIRK(a,b,c,f,x,y,x0,y0,h,x1) `): end: E0:=[1]: A2:=[3/2,-1/2]: RK4:= [[0,0,0,0],[0,0,0,1/2],[0,0,0,1/2],[0,0,0,1]],[1/6,1/3,1/3,1/6],[0,1/2,1/2,1]: GL:=[[1/4,1/4-1/6*sqrt(3)],[1/4+1/6*sqrt(3),1/4]],[1/2,1/2],[1/2-sqrt(3)/6,1/2+sqrt(3)/6]: #AAB(f,x,y,x0,y0,h,x1,L), inputs #(i) an expression f in variables x, and y #(describing the function (x,y)->f(x,y) #(ii)initial conditions y(x0)=y0 #(iii) h: mesh size #(iv): x1 last point needed #(v): L a list of numbers describing the Adams-Bashford #method #y(x+h) ~ y(x)+ h*add(L[i+1]*f(x-h*i,y[x-h*i]),i=0..nops(L)-1) #e.g. Euler is L=[1] #2-step Adams-B: L:=[3/2,-1/2] AAB:=proc(f,x,y,x0,y0,h,x1,L) local Lx,Ly,i,xc,yc: Lx:=[x0]: Ly:=[y0]: xc:=x0: yc:=y0: for i from 1 to nops(L)-1 do yc:=yc+h*subs({x=xc,y=yc},f): xc:=xc+h: Lx:=[op(Lx),xc]: Ly:=[op(Ly),yc]: od: while xcs or nops(c)<>s then RETURN(FAIL): fi: if add(b[i],i=1..s)<>1 then RETURN(FAIL): fi: if c[1]<>0 then RETURN(FAIL): fi: Lx:=[]: Ly:=[]: xc:=x0: yc:=y0: while xc<=x1 do Lx:=[op(Lx),xc]: Ly:=[op(Ly),yc]: yc:=GRK1(a,b,c,f,x,y,xc,yc,h): xc:=xc+h: od: Lx,Ly: end: ####end from C14.tx #GIRK1(a,b,c,f,x,y,x0,y0,h): Implicit #Runge-Kutte #General Runge-Kutta with Butcher matrix a,b,c #(see wiki) GIRK1:=proc(a,b,c,f,x,y,x0,y0,h) local k,s,i,eq,var,j: s:=nops(b): var:={seq(k[i],i=1..s)}: eq:= { seq(k[i]-h*subs({x=x0+c[i]*h,y=y0+add(a[i][j]*k[j],j=1..s)},f), i=1..s)}: var:=fsolve(eq,var): y0+add(b[i]*subs(var,k[i]),i=1..s): end: #GIRK(a,b,c,f,x,y,x0,y0,h,x1): The general Runge-Kutta #with Butcher matrix a,b,c appx. to #the ode #y'(x)=f(x,y(x)), y(x0)=y0 all the way #to x1 with step-size h #the output is pair of lists #Xvalues, Yavlues GIRK:=proc(a,b,c,f,x,y,x0,y0,h,x1) local xc,yc,Lx,Ly,s,i: if not (type(a,list) and type(b,list) and type(c,list)) then RETURN(FAIL): fi: s:=nops(b): if nops(a)<>s or nops(c)<>s then RETURN(FAIL): fi: if add(b[i],i=1..s)<>1 then RETURN(FAIL): fi: Lx:=[]: Ly:=[]: xc:=x0: yc:=y0: while xc<=x1 do Lx:=[op(Lx),xc]: Ly:=[op(Ly),yc]: yc:=GIRK1(a,b,c,f,x,y,xc,yc,h): xc:=xc+h: od: Lx,Ly: end: ############################################ Class Notes March 3rd ################################################ #The cheating way (using Maple) to solve the diff.eq #y'(x)=f(x,y(x)) y(x0)=y0 DS:=proc(f,x,y,x0,y0) local s: s:=dsolve({diff(y(x),x)=subs(y=y(x),f), y(x0)=y0},y(x)): if s=NULL then RETURN(FAIL): else RETURN(op(2,s)): fi: end ############################################ End of Class Notes ################################################### ############## # Problem 1a # ############## #BestAdams(f,x,y,x1) inputs an expression f in x, and y such that #dsolve({diff(y(x),x)=f,y(0)=1},y(x)); gives you a closed form answer #and outputs the d in {1,2,..,10} and h in {1/10,1/100,1/1000} such that #AAB(f,x,y,0,1,h,x1,FAB(d)) is as close as possible to the true y(x1), #and also returns the error. BestAdams:=proc(f,x,y,x1) local error1, error2, Dsolve, d, h, mesh, dstep: error1:=infinity: Dsolve:=subs(x=x1,DS(f,x,y,0,1)): for d in [1,2,3,4,5,6,7,8,9,10] do for h in [1/10, 1/100, 1/1000] do error2:=abs(Dsolve-AAB(f,x,y,0,1,h,x1,FAB(d))[2][-1]): if error2