我想寫一個使用楓樹梯形法則的程序

Question

我想寫一個使用楓樹的梯形法則程序，我會發布我的代碼，但問題是當我讀這篇文章時，我發現很難編程所有數值例子，任何人都可以幫助我嗎？

文章在這個網址：http://www.cs.technion.ac.il/~asidi/Sidi_Journal_Papers/P097_JCAM.sinpq.trans.pdf 我楓木代碼：

with(linalg); with(LinearAlgebra); with(plots); with(MTM); 
a := 0; b := 1; p := 20; 
k := 1; R := vector(p); 
f := proc (x) options operator, arrow; x*(1-x)/(1+x) end proc; 
Uex := int(f(x), x = a .. b); 
Uex := evalf(Uex); 
for N to p do 
    h := (b-a)/N; 
    U[N] := (f(a)+f(b))*(1/2); 
    for j to N-1 do 
     U[N] := U[N]+f(a+j*(b-a)/N) 
    end do; 
    U[N] := h.U[N] 
end do; 
for N to p do 
    R[N] := -evalf(log(abs(Uex-U[N]))/log(10)); 
    print(N, evalf(Uex-U[N])) 
end do; 
plot(f); listplot(R) 

Answer 1

試試這個：

a:= 0: b:= 1: p:= 20: 
R:= Vector(p): H:= Vector(p): 
f:= x-> x*(1-x)/(1+x); 
Uex:= int(f(x), x= a..b); 
for n to p do 
    N:= 2^n; 
    h:= (b-a)/N; 
    H[n]:= h; 
    R[n]:= evalhf(abs(Uex-h*((f(a)+f(b))/2 + add(f(a+j*h), j= 1..N-1))))  
end do: 
unassign('h'); 
err = Statistics:-PowerFit(H,R,h);