#Nathan Fox #Homework 18 #I give permission for this work to be posted online #Read procedures from class read(`C18.txt`): read(`hw17.txt`): Help:=proc(): print(` SIrecS(L,Ini,f,m,N) , dBVPIc(L,Ini,Fini,f,m,N)`): print(` P1dClever(f,x,h,mu0,mu1) , GlEclever(f,x,h,mu0,mu1) `): print(` EmpiricalStability(f,x,mu0,mu1,ListH) `): end: ##PROBLEM 1## #SIrecS(L,Ini,f,m,N): same input as Srec but in addition an #expression, f, in the discerte variable m, that solves #Inhomogeneous linear recurrences #a(n)=L[1]*a(n-1)+L[2]*a(n-2)+ ... + L[-1]*a(n-nops(L))+f(n) #outputs the LIST of length N+1 such that #L[i+1]=SIrecS(L,Ini,f,m,i) SIrecS:=proc(L, Ini, f, m, N) local n: return [seq(SIrec(L, Ini, f, m, n), n=0..N)]: end: ##PROBLEM 2## #dBVPIc(L,Ini,Fini,f,m,N): a clever way #to do dIBVP(...) #Discrete boundary value problem #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))+f(n) #a(0)=Ini[1],...,a(nops(Ini)-1)=Ini[-1], #a(N)=Fini[1],a(N-1)=Fini[2],... dBVPIc:=proc(L,Ini,Fini,f,m,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:=SIrecS(L, Ini1, f, m, N): eq:={seq(Li[N+2-i]=Fini[i], i=1..nops(Fini))}: var:=solve(eq, var): if var = NULL then return FAIL: else return expand(subs(var, Li)): fi: end: ##PROBLEM 3## #P1dClever(f,x,h,mu0,mu1): approximate solution of the 1-D Poisson equation #u''(x)=f(x), u(0)=mu0, u(1)=mu1 #where f is an expression in the variable x #h is the mesh size, and #h^2*u''(x) ~ u(x+h)-2*u(x)+u(x-h) #This is faster than P1d but gives the same output P1dClever:=proc(f, x, h, mu0, mu1) local m: return dBVPIc([2, -1], [mu0], [mu1], h^2*subs(x=(m-1)*h, f), m, 1/h): end: ##PROBLEM 4## #Global error in using a mesh size h #to solve the 1D Poisson equation #u''(x)=f(x), u(0)=mu0, u(1)=mu1 #Uses P1dClever instead of P1d GlEclever:=proc(f,x,h,mu0,mu1) local R1, A1, i: R1:=P1ch(f,x,h,mu0,mu1): A1:=P1dClever(f,x,h,mu0,mu1): return max(seq(abs(R1[i]-A1[i]), i=1..nops(R1))): end: ##PROBLEM 5## #EmpiricalStability(f,x,mu0,mu1,ListH): finds the equence of ratios #GlEclever(f,x,h,mu0,mu1)/h^2 #for all h in ListH EmpiricalStability:=proc(f,x,mu0,mu1,ListH) local h: return [seq(GlEclever(f,x,h,mu0,mu1)/h^2, h in ListH)]: end: #EmpiricalStability(x^4,x,0,1,[1/10,1/20,1/30,1/40]); #returned [977/25000, 188461/4800000, 29887/759375, 8061/204800]. #Using evalf yields [0.03908000000, 0.03926270833, 0.03935736626, 0.03936035156] #EmpiricalStability(exp(x),x,0,1,[1/10,1/20,1/30,1/40]); #returned [0.017526126, 0.017647884, 0.017649438, 0.01764843] (when using evalf) #EmpiricalStability(1/(3-x),x,0,1,[1/10,1/20,1/30,1/40]); took too long