Gauss Integration
>function sp(f,g) &&= integrate(f*g,x,-1,1)
integrate(f g, x, - 1, 1)
>p0 &= 1;
>p1 &= x;
>p2 &= x^2-sp(x^2,p0)/sp(p0,p0)*p0
2 1
x - -
3
>p3 &= x^3-sp(x^3,p1)/sp(p1,p1)*p1
3 3 x
x - ---
5
>p4 &= x^4-sp(x^4,p2)/sp(p2,p2)*p2-sp(x^4,p0)/sp(p0,p0)*p0 | ratsimp
4 2
35 x - 30 x + 3
-----------------
35
>fracprint(legendre(4)*8)
[ 3 0 -30 0 35 ]
>&solve(p4), xnode=%()
sqrt(2 sqrt(30) + 15) sqrt(2 sqrt(30) + 15)
[x = - ---------------------, x = ---------------------,
sqrt(35) sqrt(35)
sqrt(15 - 2 sqrt(30)) sqrt(15 - 2 sqrt(30))
x = - ---------------------, x = ---------------------]
sqrt(35) sqrt(35)
[ -0.861136311594 0.861136311594 -0.339981043585 0.339981043585 ]


>A=xnode^((0:3)');
>b=(2*[1,0,1/3,0])';
>alpha=(A\b)'
[ 0.347854845137 0.347854845137 0.652145154863 0.652145154863 ]
>function map gauss5test (f) := sum(f(xnode)*alpha)
>fracprint(gauss5test(["1","x","x^2","x^3","x^4","x^5","x^6","x^7","x^8"]))
[ 2 0 2/3 0 2/5 0 2/7 0 258/1225 ]
>function gauss5test (f,a:number,b:number,n=1,xn=xnode,al=alpha) ...
h=(b-a)/n;
if n>1 then x=linspace(a+h/2,b-h/2,n-1); else x=a+h/2; endif;
X=x-(h/2*xn)';
return sum(al.f(X))*h/2
endfunction
>fracprint(gauss5test("x^6",0,1))
1/7
>longestformat; gauss5test("exp(x)",0,1,10)-(E-1)
4.440892098500626e-016
>longestformat; gauss5("exp(x)",0,1,10)-(E-1)
4.440892098500626e-016

>simpson("exp(x)",0,1,10)-(E-1)
5.964481175624314e-008
>gauss("exp(x)",0,1)-(E-1)
2.220446049250313e-016
>gauss5("exp(-x^2/2)",0,10,n=10)/sqrt(2*pi)
0.5000000000028569
>gauss("exp(-x^2/2)",0,10,n=10)/sqrt(2*pi)
0.5000000000000001
>function sp(f,g) &&= integrate(f*g*exp(-x),x,-1,1)
- x
integrate(f g E , x, - 1, 1)
>p0 &= 1;
>p1 &= x-sp(x,p0)/sp(p0,p0)|ratsimp
2
(E - 1) x + 2
--------------
2
E - 1
>p2 &= x^2-sp(x^2,p1)/sp(p1,p1)*p1-sp(x^2,p0)/sp(p0,p0)*p0|ratsimp
4 2 2 4 2 4 2
(E - 6 E + 1) x + (- 2 E + 16 E - 6) x - E + 6 E + 7
-----------------------------------------------------------
4 2
E - 6 E + 1
>function p(x,i) := x^i*exp(-x)
>function map a(i) := gauss("p",-1,1,100;i)
>a(5)
-0.3242973696922066
>&integrate(x^5*exp(-x),x,-1,1), %()
- 1
44 E - 326 E
-0.3242973696922178
>A=a((0:10)+(0:9)')_1; b=zeros(10,1)_1;
>p=xlgs(A,b)';
>plot2d("polyval(p,x)",-1,1); insimg;

>for i=0 to 9; gauss("x^i*polyval(p,x)*exp(-x)",-1,1,100), end;
1.23084875625068e-014
-7.745415420146173e-014
-1.233568802661011e-014
3.418543226274551e-014
3.103073353827313e-016
-2.974398505273257e-014
-1.885935851930754e-014
6.401379426534959e-014
3.972544515562504e-014
-5.778488798569015e-014
>xn=sort(real(polysolve(p)))
Column 1 to 2:
-0.9761983950159177 -0.876341759205208
Column 3 to 4:
-0.7038368925296208 -0.4708613994245012
Column 5 to 6:
-0.1948808338591103 0.1018814014714706
Column 7 to 8:
0.3935411401555917 0.6524949822560991
Column 9 to 10:
0.8523013113814679 0.9712722677116907
>M=xn^((0:9)'); b=a((0:9)');
>al=xlgs(M,b)'
Column 1 to 2:
0.1616042035976784 0.330895277283453
Column 3 to 4:
0.414855139626519 0.4128502394503029
Column 5 to 6:
0.3529200553213236 0.269811110073977
Column 7 to 8:
0.1888527845927686 0.1216033541103824
Column 9 to 10:
0.06925902427714206 0.02775119895405546
>function gauss10exp (f,a:number,b:number,n=1,xn=xn,al=al) ...
h=(b-a)/n;
if n>1 then x=linspace(a+h/2,b-h/2,n-1); else x=a+h/2; endif;
X=x-(h/2*xn)';
return sum(al.f(X))*h/2
endfunction
>gauss10exp("cos(x)",-1,1)
1.933421496200713
>sum(cos(xn)*al)
1.933421496200713
>&integrate(cos(x)*exp(-x),x,-1,1); %()
1.933421496200713
Examples Homepage