#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: