#hw23.txt, 19 April 2014, Andrew Lohr Pw:=proc(Lf,x,Lx) local i,command: command:= x=Lx[i-1] and x=Lx[nops(Lx)],0): end: Compare :=proc(p,q,v,x,x1) local u,f,N,ret: u := v*x*(1-x): f:=-diff(p*diff(u,x),x) + q*u: ret:=[]: for N in {10,20,30,40} do ret:=[op(ret),RR(p,q,f,x,[seq(i/N,i=0..N)])[round(x1*N+1)]]: od: evalf(subs(x=x1,[op(ret),u])): end: #running it with # Compare(x^2,1+x^3,sin(x),x,.513); # Compare(x^2,1+x^3,x^7+1,x,.513); # Compare(x^2+3,sin(x),exp(x),x,.542); #we get # [0.1251017497, 0.1234668473, 0.1231112603, 0.1227372190, 0.1226153759] # [0.2562016786, 0.2535615542, 0.2529879707, 0.2522113970, 0.2521669733] # [0.4355599517, 0.4264488333, 0.4274669885, 0.4266447327, 0.4268274813] BSpline:=proc(d,Lx,i) local a,Lu,i1,j1,var,eq: if nops(Lx)<>d+2 then RETURN(FAIL): fi: Lu:=[seq(add(a[i1,j1]*x^j1,j1=0..d),i1=0..d)]: var:={seq(seq(a[i1,j1],j1=0..d),i1=0..d)}: eq:={subs(x=Lx[1],Lu[1]), seq(subs(x=Lx[1],diff(Lu[1],x\$i1)),i1=1..d-1)}: eq:= eq union {subs(x=Lx[d+2],Lu[d+1]), seq(subs(x=Lx[d+2],diff(Lu[d+1],x\$i1)),i1=1..d-1)}: #must be 1 at proscribed position eq:= eq union {subs(x=Lx[i+1],Lu[i] = 1)}: #smoothness conditions eq:= eq union {seq(subs(x=Lx[i1],Lu[i1] = Lu[i1-1]),i1 = 2..d+1)}: eq := eq union {seq( seq( subs(x=Lx[i1],diff(Lu[i1-1],x\$j1)) - subs(x=Lx[i1],diff(Lu[i1],x\$j1)),j1=1..d-1) ,i1= 2..d+1)}: var:=solve(eq,var): if(var = NULL) then return(FAIL): fi: Pw(subs(var,Lu),x,Lx): end: TestF3:=proc(c,Lx,x) local i: expand(add(c[i]*BSpline(nops(Lx)-2,Lx,i+1),i=1..nops(Lx)-2)), { seq( c[i], i=1..nops(Lx)-2)}: end: RR3:=proc(p,q,f,x,Lx) local c, Lu, F, eq,var: Lu:=TestF3(c,Lx,x): var:=Lu[2]: Lu:=Lu[1]: F:=EvalFunct(p,q,f,x,Lx,Lu): eq:=evalf({seq(diff(F,var[i]),i=1..nops(var))}): var:=solve(eq,var): if var=NULL then FAIL: else subs(var,Lu): fi: end: