#!/usr/local/bin/maple # -*- maplev -*- # Nathaniel Shar # HW 17 - 30 Mar 2014 # Experimental Mathematics ############# # Problem 1 # ############# SIrec := proc(L, initial, f, N) local i: option remember: if N < nops(initial) then: return initial[N+1]: else: return add(L[i]*SIrec(L, initial, f, N-i), i=1..nops(L)) + f(N): fi: end: SIrecS := proc(L, initial, f, N) local i: return [seq(SIrec(L, initial, f, i), i=0..N)]: end: ############# # Problem 2 # ############# dBVPIc := proc(L, initial, final, f, N) local eqns, vars, solns, A, newinitial, i, R: if nops(L) <> nops(initial) + nops(final) then: return FAIL: fi: vars := {seq(A[i], i=1..nops(final))}: newinitial := [op(initial), op(vars)]: R := SIrecS(L, newinitial, f, N): eqns := {seq(R[-i] = final[i], i=1..nops(final))}: solns := solve(eqns, vars): if solns = NULL then: return FAIL: fi: return subs(solns, R): end: ############# # Problem 3 # ############# P1dClever := proc(f, h, mu0, mu1) local N: N := floor(1/h): return dBVPIc([2, -1], [mu0], [mu1], x->h^2*f(h*(x-1)), N): end: ############# # Problem 4 # ############# 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: GIEClever := proc(f, h, mu0, mu1) local A, B, x: A := P1ch(f(x), x, h, mu0, mu1): B := P1dClever(f, h, mu0, mu1): return max(seq(abs(A[i] - B[i]), i=1..nops(B))): end: ############# # Problem 5 # ############# EmpiricalStability := proc(f, mu0, mu1, L): return [seq(GIEClever(f, h, mu0, mu1)/h^2, h=L)]: end: # x^4: 0.03908000000, 0.03926270833, 0.03935736626, 0.03936035156 # exp: 0.017526126, 0.017647883, 0.017649542, 0.017649353 # 1/(3-x): 0.001387500000, 0.001402700000, 0.001401000000, 0.001403000000