kill(all)$

method :undetermined coefficients, 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(f(x),x,a,b) = sum(c[i]*f(x[i]),i,0,n);

a) suppose : a ≤ x[i] ≤ b
b) c[i]:undetermined coefficients

example : a = 0, b=h ,x[0]=0,x[1]=h/3,x[2]=h

1e) step : given

p2a:[a=0,b=h,n=2];

→ domain : x[0] ≤ x[1] ≤ x[2]  → 0≤h/3≤h,unequal spaced

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,2);

→values

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

p2b1a1:[0,h/3,h];

p2b1a2:map("=",p2b1a,p2b1a1);

3e) step:find equations from the polynomials.

need  polynomial:f(x) ={1,x,x^2,..,x^(n-1)}

p2c:makelist(x^i,i,0,2);

p2c1:makelist(f(x)=p2c[i],i,1,3);

lhs :0,h/3,h

powerdisp:true;

→ first equation

p2d:makelist(f(x[i])=(x[i])^0,i,0,2);

p2d1:lratsubst(p2d,rhs(p2b))=subst([p2c1[1]],lhs(p2b));

→ all equations

p2e:s[k]:=makelist(f(x[i])=(x[i])^k,i,0,2);

 => first equation

s[0];/*same p2d */

lratsubst(s[0],rhs(p2b))=subst([p2c1[1]],lhs(p2b));

info : s[0]  →  p2c1[1]

 => second equation

s[1];

p2f1:lratsubst(s[1],rhs(p2b))=subst([p2c1[2]],lhs(p2b));

info : s[1]  →  p2c1[2]

p2f2:ev(lratsubst(p2b1a2,p2f1),integrate);

 =>  all

p2g[k]:=lratsubst(s[k],rhs(p2b))=subst([p2c1[k+1]],lhs(p2b));

info : s[k]  →  p2c1[k+1]   , k =0,1,..,n

p2ga:makelist(p2g[k],k,0,2);

p2gb:ev(lratsubst(p2b1a2,p2ga),integrate);

4e) step:find c[i], i=0,1,..n.

p2h:solve(p2gb,p2b11a);

p2ha:flatten(p2h);

5e) step:find integration equation

p2i:p2ha;

p2i1:p2b1a2;

p2i2:flatten(cons(p2i,p2i1));

p2i3:factor(lratsubst(p2i2,p2b));

find :error ?

try :1e) f(x): =x^2

p3a:f(x):=x^2;

p3b:error1 = ev(lhs(p2i3)-rhs(p2i3),nouns);

forget(f); /*forget f(x) */

p4a:f(x):=x^3;

p3b:error2 = ev(lhs(p2i3)-rhs(p2i3),nouns);

method exact :polynomials :{x^0,x^1,x^2}

kill(f);

integrate_error = rhs(p3b)*at('diff(f(x),x,3),x=ζ);

info : domain :0< ζ < h