#hw16.txt Ross Berkowitz # read("E:\\Documents\\Maple\\em14\\hw16.txt"); ############# 1 ################# #BestAdams(f,x,y,x1) inputs an expression f in the variables x and y such that #dsolve({diff(y(x),x)=f,y(0)=1},y(x)) works. Then tells you for which choices of #d between 1 and 10 and h in {1/10,1/100,1/1000} such that AAB(f,x,y,0,1,x1,FAB(d)) #is as close to the truth as possible and returns the error BestAdams:=proc(f,x,y,x1) local truth,hopeful,hope,champ,record,d,h,D,H,ans: ans:=rhs(dsolve({diff(y(x),x)=subs(y=y(x),f), y(0)=1},y(x))): truth:=evalf(subs(x=x1,ans)): if ans= null then print(`Expression is two complicated to solve, try again`): RETURN(FAIL): fi: record:=infinity: H:={1/10,1/100,1/100}: D:={seq(d,d=1..10)}: for h in H do for d in D do hope:=abs(AAB(f,x,y,0,x1,h,x1,FAB(d))[2][-1]-truth): if hopef(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: