#March 31, 2014, C18.txt difference methods Help:=proc(): print(` Srec(L, Ini, n), dBVP(L,Ini,Fini,N) `): print(`SrecS(L,Ini,N) , dBVPc(L,Ini,Fini,N), P1ch(f,x,h,mu0,mu1)`): print(`P1d(f,x,h,mu0,mu1), GlE(f,x,h,mu0,mu1) `): end: #Srec(L,Ini,n): The n-th term of the #solution of the linear recurrence equation #with constant coefficients #a(n)=L[1]*a(n-1)+L[2]*a(n-2)+ ... + L[-1]*a(n-nops(L)) #a(0)=Ini[1], ..., a(nops(Ini)-1)=Ini[-1] Srec:=proc(L,Ini,n) local i: option remember: if nops(L)<>nops(Ini) then print(Ini, L, `should be lists of the SAME length `): RETURN(FAIL): fi: if n<0 then 0: elif nnops(Ini)+nops(Fini) then RETURN(FAIL): fi: var:={seq(a[i],i=0..N)}: eq:={seq(a[i]=Ini[i+1],i=0..nops(Ini)-1), seq(a[N-i]=Fini[i+1],i=0..nops(Fini)-1)}: eq:= eq union {seq(a[i]-add(L[j]*a[i-j],j=1..nops(L)), i=nops(L)..N)}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): else subs(var,[seq(a[i],i=0..N)]): fi: end: #####end old stuff #SrecS(L,Ini,N): inputs a linear recurrence equation #with constant coefficients, coded by L #and a list of initial values, Ini #outputs the LIST of length N+1 such that #L[i+1]=Srec(L,Ini,n) SrecS:=proc(L,Ini,N) local n: [ seq(Srec(L,Ini,n),n=0..N)]: end: #dBVPc(L,Ini,Fini,N): a (hopefully) clever way #to do dBVP(...) #The list of length N+1 #giving the values of a(n) for n=0 to n=N #to the solution of the linear recurrence equation #with constant coefficients #a(n)=L[1]*a(n-1)+L[2]*a(n-2)+ ... + L[-1]*a(n-nops(L)) #a(0)=Ini[1], ..., a(nops(Ini)-1)=Ini[-1]. #a(N)=Fini[1], a(N-1)=Fini[2], .. dBVPc:=proc(L,Ini,Fini,N) local eq,var,a,i,Ini1, Li: Ini1:=[op(Ini),seq(a[i],i=1..nops(Fini))]: var:={seq(a[i],i=1..nops(Fini))}: Li:=SrecS(L,Ini1,N): #Li[N+1]=Fini[1], Li[N]=Fini[2], Li[N-1]=Fini[3]: eq:={seq(Li[N+2-i]=Fini[i], i=1..nops(Fini))}: var:=solve(eq,var): subs(var,Li): end: #P1ch(f,x,h,mu0,mu1): the exact solution of the 1-D Poisson eq. #u''(x)=f(x), u(0)=mu0, u(1)=mu1 #where f is an expression in the variable x P1ch:=proc(f,x,h,mu0,mu1) local u,c1,c2,eq,var: u:=int(int(f,x)+c1,x)+c2: var:={c1,c2}: eq:={subs(x=0,u)=mu0, subs(x=1,u)=mu1}: var:=solve(eq,var): u:=subs(var,u): [seq(subs(x=i*h,u),i=0..1/h)]: end: #P1d(f,x,h,mu0,mu1): approximate solution of the 1-D Poisson eq. #u''(x)=f(x), u(0)=mu0, u(1)=mu1 #where f is an expression in the variable x #h the mesh size, and #u''(x) ~ u(x+h)-2*u(x)+u(x-h) P1d:=proc(f,x,h,mu0,mu1) local u,Li,var,eq: Li:=[seq(u[i],i=0..1/h)]: var:={seq(u[i],i=0..1/h)}: eq:={ seq(u[i+1]-2*u[i]+u[i-1]=h^2*subs(x=i*h,f), i=1..1/h-1)}: eq:=eq union {u[0]=mu0, u[1/h]=mu1}: var:=solve(eq,var): subs(var,Li): end: #Global error in using a mesh size h #to solve the 1D Poisson eq. #u''(x)=f(x), u(0)=mu0, u(1)=mu1 GlE:=proc(f,x,h,mu0,mu1) local R1, A1: R1:=P1ch(f,x,h,mu0,mu1): A1:=P1d(f,x,h,mu0,mu1): max(seq(abs(R1[i]-A1[i]),i=1..nops(R1))): end: