pdf_numerik

179 days ago by PatrickHammer

def gauss_legendre(f,n,a,b): P=lambda n:1/(2^n*factorial(n))*diff(((x^2-1)^n),x,n) X=[u.rhs() for u in solve(P(n+1)==0,x)] scale(t)=(t*(b-a)+a+b)/2 return (b-a)/2*sum([integrate(prod([(x-X[j])/(X[i]-X[j]) for j in [0,..,n] if i!=j]),x,-1,1)*f(scale(X[i])) for i in [0,..,n]]) def newton_cotes(f,n,a,b): X=lambda i:a+(b-a)/n*i return sum([integrate(prod([(x-X(j))/(X(i)-X(j)) for j in [0,..,n] if i!=j]),x,a,b)*f(X(i)) for i in [0,..,n]]) f(x)=sin(x); show(f) l=gauss_legendre(f,4,-1,2).n() m=newton_cotes(f,4,-1,2).n() e=f.integrate(x,-1,2).n() print "legendre:",l,"milne:",m,"exact:",e,"error_legendre:",abs(l-e),"error_milne:",abs(m-e) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ \sin\left(x\right)

legendre: 0.956449174395270 milne: 0.955949089353632 exact:
0.956449142415282 error_legendre: 3.19799882131377e-8 error_milne:
0.000500053061650130
\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ \sin\left(x\right)

legendre: 0.956449174395270 milne: 0.955949089353632 exact: 0.956449142415282 error_legendre: 3.19799882131377e-8 error_milne: 0.000500053061650130
def trapez_sum(f,n,a,b): h=(b-a)/n return h*((f(a)+f(b))/2+sum([f(a+c*h) for c in range(1,n)])) f=1/sqrt(1-1/2*sin(x-1)); show(f) for a in [0,1]: for N in [2^i for i in [0,..,5]]: v=n(trapez_sum(f,N,a,2*pi)) w=n(f.integrate(x,a,2*pi)) print "a",a,"N",N,"exact:",w,"trapez:",v,"error:",abs(v-w) 
        
\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{\sqrt{-\frac{1}{2} \, \sin\left(x - 1\right) + 1}}

a 0 N 1 exact: 6.62655268095 trapez: 5.27136699864820 error:
1.35518568229818
a 0 N 2 exact: 6.62655268095 trapez: 6.76341474604790 error:
0.136862065101520
a 0 N 4 exact: 6.62655268095 trapez: 6.61414918603767 error:
0.0124034949087077
a 0 N 8 exact: 6.62655268095 trapez: 6.62654244910166 error:
0.0000102318447199323
a 0 N 16 exact: 6.62655268095 trapez: 6.62655267967001 error:
1.27636745617110e-9
a 0 N 32 exact: 6.62655268095 trapez: 6.62655268094638 error:
1.77635683940025e-15
a 1 N 1 exact: 5.72129699514 trapez: 4.85779409557072 error:
0.863502899565689
a 1 N 2 exact: 5.72129699514 trapez: 5.45843920271536 error:
0.262857792421049
a 1 N 4 exact: 5.72129699514 trapez: 5.69888127071675 error:
0.0224157244196528
a 1 N 8 exact: 5.72129699514 trapez: 5.71509258986420 error:
0.00620440527220367
a 1 N 16 exact: 5.72129699514 trapez: 5.71974941218166 error:
0.00154758295474622
a 1 N 32 exact: 5.72129699514 trapez: 5.72091025218410 error:
0.000386742952306562
\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{\sqrt{-\frac{1}{2} \, \sin\left(x - 1\right) + 1}}

a 0 N 1 exact: 6.62655268095 trapez: 5.27136699864820 error: 1.35518568229818
a 0 N 2 exact: 6.62655268095 trapez: 6.76341474604790 error: 0.136862065101520
a 0 N 4 exact: 6.62655268095 trapez: 6.61414918603767 error: 0.0124034949087077
a 0 N 8 exact: 6.62655268095 trapez: 6.62654244910166 error: 0.0000102318447199323
a 0 N 16 exact: 6.62655268095 trapez: 6.62655267967001 error: 1.27636745617110e-9
a 0 N 32 exact: 6.62655268095 trapez: 6.62655268094638 error: 1.77635683940025e-15
a 1 N 1 exact: 5.72129699514 trapez: 4.85779409557072 error: 0.863502899565689
a 1 N 2 exact: 5.72129699514 trapez: 5.45843920271536 error: 0.262857792421049
a 1 N 4 exact: 5.72129699514 trapez: 5.69888127071675 error: 0.0224157244196528
a 1 N 8 exact: 5.72129699514 trapez: 5.71509258986420 error: 0.00620440527220367
a 1 N 16 exact: 5.72129699514 trapez: 5.71974941218166 error: 0.00154758295474622
a 1 N 32 exact: 5.72129699514 trapez: 5.72091025218410 error: 0.000386742952306562