kill(all)$

Maxima :Gaussian Quadrature ,version 1.0
peter.vlasschaert@gmail.com,02/10/2019

load(lrats)$

How to find approximation : d1

d1:integrate(f(t),t,-1,1);

1e)  Gauss-Legendre (Two-point rule)
=============================

eq1:a1+a2=2;                        /*f(x) =1 */
eq2:a1*x1+a2*x2=0;             /*f(x) = x */
eq3:a1*x1^2+a2*x2^2=2/3; /*f(x) = x^2 */
eq4:a1*x1^3+a2*x2^3=0;    /*f(x) = x^3 */

rhs: {eq1,eq2,eq3,eq4}

q1:makelist(integrate(t^i,t,-1,1),i,0,3);

solve([eq1,eq2,eq3,eq4],[a1,a2,x1,x2]);

2e)  Gauss-Legendre (Three-point rule)
==============================

eq1a:a1+a2+a3=2;
eq2a:a1*x1+a2*x2+a3*x3=0;
eq3a:a1*x1^2+a2*x2^2+a3*x3^2=2/3;
eq4a:a1*x1^3+a2*x2^3+a3*x3^3=0;
eq5a:a1*x1^4+a2*x2^4+a3*x3^4=2/5;
eq6a: x3=0;                       /*missing equation */

rhs: {eq1a,eq2a,eq3a,eq4a,eq5a}

q2:makelist(integrate(t^i,t,-1,1),i,0,4);

s1:solve([eq1a,eq2a,eq3a,eq4a,eq5a,eq6a],[a1,a2,a3,x1,x2,x3]);

s1a:s1[1];
s1b:s1[2];

p1:f(x):=1/(x+3);

How make calculation:

p1i : 'integrate(f(t),t,-1,1);

f(0) = ?

p1a:lratsubst([x1=0],f(x1));

p2:lratsubst(s1a,a1*f(x1)+a2*f(x2)+a3*f(x3));

approximation:integral , p1i
=====================

p2a:float(p2);

exact result
=========

p2b:p1i=ev(p1i,nouns);

p2c:float(p2b);

error : value 'p1i'
=========

p2d:error = rhs(p2c) -p2a;

plot : f(x)

[parametric, r*sin(t), r*cos(t), [t, -8*%pi, 8*%pi]

wxplot2d([f(x),[parametric, 1,x, [x, 0, f(1)]]], [x,-2,2],[y,0,1]);

t1:1;
t2:-1;

see :website:grahics ,(truebasic,maxima)

w1a:f(x);                              /*function*/
w1b:[parametric, t1,x, [x, 0, f(t1)]];  /*vertical lines crossing the function */
w1c:[parametric, t2,x, [x, 0, f(t2)]];

inline plot2d

wxplot2d([w1a,w1b,w1c ],[x,-1+t2,+1+t1],[y,0,1],[color,blue,red,red]);

rem :value of the integral is the area below the function into interval .
===================================================
x: t2 ≤ x ≤ t1
f(x): 1/(x+3)

p3:'integrate(f(t),t,-1,1) =integrate(f(t),t,-1,1) ;

p3a:print( "area of the integral f(x) = ",p2a,":","range =>",t2," ≤  x ≤",t1,"")$

?  What to do : integration range different from above
range : a ≤y ≤ b  , a,b are real values

1e) solution : linear transformation 
==========================

p4:x=(b-a)/2*y+(b+a)/2;

2e) differentation : linear transformation 
==============================

declare(a,constant,b,constant);

p4a:diff(p4);

p4a1:lhs(p4a);
p4a2:rhs(p4a);
p4a3:coeff(p4a2,del(y));

3e) conversion
============
general : f(x) = g(x)

p5:integrate(g(x),x,a,b) = p4a3*integrate(g(rhs(p4)),y,-1,1);

4e) general gaussian quadrature (= approximation)
=======================================

p5a:lhs(p5)= p4a3*sum(A[i]*g((b-a)/2*x[i]+((b+a)/2)),i,1,n);

rem : 'y range' => ' x[i] range '

we start discussion : a1,a2,a3,  .....  => A[i]
                                    x1,x2;x3, .........=> x[i]
n=2: Gauss-Legendre : Two-point rule
n=3 :Gauss-Legendre:  Three-point rule

same computation : general  equation , 'p5a' for f(x) integration
from above. 
================================================

assume(b>a);

1e) step : g => f
-----------------------

p6a:subst([g=f],p5a);

2e) values a,b
----------------------

p6b:[a=-1,b=1];

p6b1:lratsubst(p6b,p6a);

3e) example: integration of f(x) =1/(x+3), -1≤ x ≤1
for n:3, three-point rule
-------------------------------------------------------------------

a)  set n=3
---------------

p6c:subst([n=3],p6b1);

p6c1:rhs(p6c),simpsum;

rem:sum expand => simpsum

p6c2:lhs(p6c) = p6c1;

b) subst values {A[i],x[i]} i=1,2,3
-------------------------------------------
use : 'above',s1a and p6c2

p6d1:[a1=A[1],a2=A[2],a3=A[3],x1=x[1],x2=x[2],x3=x[3]];
p6d2:lratsubst(p6d1,s1a);

p6d3:lhs(p6c2)=lratsubst(p6d2,rhs(p6c2));

p6d4:float(p6d3);

 =========================================================