kill(all)$

maxima:gauss-legendre quadrature, integration,version 1.0
peter.vlasschaert@gmail.com, 23/12/2020

load packages,option
----------------------------

load(lrats)$ ratprint:false $

concept: integration , 'finite domain', [a,b] ,@ a < b
=======================================
given : f(x) = polynomial

 approximation of the integral : interval  [a,b]
------------------------------------------------------------

declare(n,integer)$

assume(n>0,a<b);

p1:integrate(w(x)*f(x),x,a,b) = sum(c[i]*f(x[i]),i,1,n);

a) suppose : a ≤ x[i] ≤ b
b) c[i]:undetermined coefficients : n
c) x[i]:undetermined                    : n
d) w(x):weigthing function 
------------------------------------------------------
exact :polynomials   : 2*n -1  → x^(2*n-1)

two-point formula: gauss-legendre quadrature

n:2;

1e) step : given

p2a:[n=2,w(x)=1];

powerdisp:true$

2e) step:approximation integral

p2b:ev(lratsubst(p2a,p1),simpsum);

p2b1:sort(listofvars(rhs(p2b)));

→ coef.

p2b11a:makelist(p2b1[2*k+1],k,0,n-1);

→ undetermined.

p2b1a:makelist(p2b1[2*k],k,1,n);

n=2,legendre polynomials
----------------------------------
legendre differential equation :
(1-x^2)*y'' -2*x*y'+n*(n+1)y =0
sol : Rodrigue's Formula

def:Rodrigue's Formula
==================

q1:P[nn](xx):=(1/(nn!*2^nn))*diff((xx^2-1)^nn,xx,nn);

q2:ratsimp(makelist('P[i](x)=P[i](x),i,0,3));

roots: P[i](x) =0 i=0,1,..,n-1
=====================
example :  n=4

→ numerical values 'roots'
------------------------------------

q3:makelist(allroots(rhs(q2[i])),i,1,length(q2));

→ exact values      'roots'
------------------------------------

→ P[0](x) =0

q3a:solve(rhs(q2[1])=0,x);

→ P[1](x) =0

q3b:solve(rhs(q2[2])=0,x);

→ P[2](x) =0

q3c:solve(rhs(q2[3])=0,x);

→ P[3](x) =0

q3d:solve(rhs(q2[4])=0,x);

Orthogonality Polynomials ( legendre (w(x)=1))
====================

orthocomb(m1,n1,ww):=integrate(ww*P[m1](x)*P[n1](x),x,-1,1);

? test   m1,n1  

first property :
===========

first     :  n1=m2=3
-------------------------
n1=m2  → orthocomb(m1,n1,ww) > 0

orthocomb(2,2,1);

second :  n1≠m2
-------------------------
n1≠m2  → orthocomb(m1,n1,ww) = 0

orthocomb(2,3,1);

? test   x^j,n1

second property :
=============
check : j < n1 => orthocomb(j,n1,ww) = 0

orthocomp(j,n1,ww):=integrate(ww*P[n1](x)*x^j,x,-1,1);

n1=3,j=2

'orthocomp(2,3,1) = orthocomp(2,3,1);

two-point formula
=============

first step    : use p2b
----------------------------

q4:p2b;

q41:subst([a=-1,b=1],q4);

second step : find list variables ,solve equations
-----------------------------------------------------------------

q4a:sort(listofvars(rhs(q4)));

thrid     step : polynomials 
----------------------------------
f(x)=[1,x,....,x^(2*n-1)]

q4b:makelist(x^i,i,0,2*n-1);

fourth  step:  2*n equations
------------------------------------

→ first equation 

q4c1:[f(x)=q4b[1],f(x[1])=q4b[1],f(x[2])=q4b[1]];

q4c1a:reverse(lratsubst(q4c1,q41));

→ second equation

q4c1b1:[f(x)=q4b[2],f(x[1])=subst([x=x[1]],q4b[2]),
    f(x[2])=subst([x=x[2]],q4b[2])];

q4c1b2:reverse(lratsubst(q4c1b1,q41));

q4c1b3:ev(%,integrate);

→ all  equation

sw[k]:=[f(x)=q4b[k],f(x[1])=subst([x=x[1]],q4b[k]),
    f(x[2])=subst([x=x[2]],q4b[k])];

sw[1];/* same as q4c1 */

qd1:eql(kk):=reverse(lratsubst(sw[kk],''q41));

qd2:makelist(eql(i),i,1,2*n);

qd3:ev(qd2,integrate);

qd3a:qd3[1];
qd3b:qd3[2];
qd3c:qd3[3];
qd3d:qd3[4];

third step : solve equations

qe1:solve(qd3,q4a);

qe1a:qe1[1];
qe1b:qe1[2];

fifth step : final result (q41)

qf1:lratsubst(qe1a,q41);
qf2:lratsubst(qe1b,q41);
qf:is(rhs(qf1)=rhs(qf2));

→ P[n](x) =0, see : n=2,q3c, 'see above'

qf3:q3c;

roots: legendre polynomials (n=2, reason : two-point method)

polynomials : legendre .
==================

representation: d(x),only symbolic (use not :dx)

qqa:integrate('P[2](x)*d(x),x,-1,1)=cc[1]*'P[2](x1)+cc[2]*'P[2](x2);

find eq1: property orthogonality

qqb1:integrate('P[2](x)*'P[0](x)*d(x),x,-1,1)=rhs(qqa);

qqb2:subst([d(x)=1],qqb1);

qqb3:ev(lhs(qqb2),nouns)=rhs(qqb2);

→ eq1

qqb4:reverse(qqb3);

from :qd3a → (1),qd3b→(2),qd3c→(3),qd3d →(4).

info: same ,'see above'
-------------------------------

find eq2: property orthogonality

ra:q2[3];

ra1:reverse(subst([x=x[1]],ra));
ra2:reverse(subst([x=x[2]],ra));

ra3:[ra1,ra2];

property : integral , odd function , symmetric (ref :origin)

--------------------------------------------------------------------------

ra1:integrate(3/2*x^3,x,-1,1);
ra2:integrate(-1/2*x,x,-1,1);

--------------------------------------------------------------------------

equations :(4),(2)

rb1:integrate(x^3,x,-1,1)=lhs(qd3d);
rb2:integrate(x,x,-1,1)=lhs(qd3b);

equations :3/2*'(4)' & -1/2*'(2)'

rb1a:3/2*rb1;
rb2a:-1/2*rb2;

rbb:expand(rb1a+rb2a);

rbb1:lhs(rbb)=collectterms(rhs(rbb),c[1],c[2]);

rbb2:lhs(rbb)=factorsum(rhs(rbb1));

rbb3a1:part(num(rhs(rbb2)),1);
rbb3a2:part(num(rhs(rbb2)),2);
rbb3b:denom(rhs(rbb2));

rbbb:lhs(rbb)=rbb3a1/rbb3b+rbb3a2/rbb3b;

→ eq2

→ use :P[2](x)=(-1+3*x^2)/2, x=x[i], i=1,2

rbbbb:lratsubst(ra3,rbbb);

result :eq1 & eq2

result : select c[i], i=1,2 arbitrary
→ P[2](x), x=x[1],x[2] be roots of polynomial

How to use : above concept

given :

first step : integral

rx1:integrate(f(Y)*d(Y),Y,a,b)= A;

second step : convert domain :[a,b] → [-1,1]

rx2:Y = (b-a)/2*t+(a+b)/2; /*linear transformation*/

check : -1,1

rx2a:expand(subst([t=-1],rx2));
rx2b:expand(subst([t=1],rx2));

rx3:'diff(Y,t)=diff(rhs(rx2),t);

rx3a:subst([a=-1,b=1],rx1);

rx4:[d(Y)=rhs(rx3),rx2];

rx5:ratsimp(lratsubst(rx4,rx3a));

A = result integral

final result :integral
===============

ry1:rhs(rx2);
ry1a:subst([t=rhs(q3c[1])],ry1);
ry1b:subst([t=rhs(q3c[2])],ry1);

ry2a:rx2;
ry2a:d(Y)=rhs(rx3);

lhs(rx1)=rhs(rx3)*(f(ry1a)+f(ry1b));