kill(all)$

maxima : simple linear integral equation ( = numerical),part 2,version 1.0
peter.vlasschaert@gmail.com, 10/03/2022

1e) ref:maxima : simple linear integral equation ( = exact),part 1,version 1.0
2e) ref:https://en.wikipedia.org/wiki/Simpson%27s_rule
see :website:maxima: appendix: smpson1/3 rule for integration,version 1.0

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

load(lrats)$ratprint:false$

given : simple integral equation
-----------------------------------------

p1a1:y(x) =1+'integrate('x*t*y(t),t,0,1);

see :1e)  exact  solution : "p1a1"

p1a2:y(x)=1+3/4*x;

solve : numerical

→ first step : divided interval :[0,1] , two parts

a:0$  /*begin of the interval: integral*/
b:1$  /*end    of  the interval: integral*/
nu:2$ /*divide number of parts*/

p2b1:[a,b]; /*interval: integral*/

p2b2:(b-a)/nu;/*step value:integral*/

values : discrete points on the interval :[a,b]

p2b:makelist(a+p2b2*k,k,0,nu);

rem : divide part of "integral part:domain" of the : Integral Equation

→second step: analyse 

p2c11:args(p1a1);

p2c12a:args(p2c11[2]);

p2c12b:args(p2c12a[1]);

p2c12c:args(p2c12b[1]); /* integrandum,var,begin:a,end:b*/

→ thrid step: discrete values, use :"p2c12c"

p2d1:makelist(subst([p2c12c[2]=p2b[k]],p2c12c[1]),k,1,length(p2b));

variables : need to be calculated , "unknown"

p2d2:map(y,p2b);

p2d3:p2c12c[2];

→ fourth step:  approximation integral

simpson rule 1/3 : h/3(1*fo,4*f1,1*f2), see 2e

integrate(f(x),x,a,b) =h/3*(f0+4*f1+2*f2+4*f3+2f4+.....+fn),here n=nu
--------------------------------------------------------------------------------
integrate(f(x),x,a,b) =h/3*(f0+4*f1+f2)
--------------------------------------------------------------------------------
rem :(1,4,1)                              → h/3*(1*f0+4*f1+1*f2)
rem :(1,4,1),(1,4,1)  = (1,4,2,4,1) → h/3*(1*f0+4*f2+2*f3+4*f4+f5)
--------------------------------------------------------------------------------

p2e1a:makelist(concat(f,i),i,0,2);

use:p2d1

p2e1b1:map("=",p2e1a,p2d1);

use : p2b2

p2e1b2:[h=p2b2];

p2e1b3:flatten([p2e1b1,p2e1b2]);

p2e1c:p2c12b[1]=expand(subst(p2e1b3,h/3*(f0+4*f1+f2)));

→fifth step : subst , "p1a1"

p3a1:subst([p2e1c],p1a1);

→sixth step :construct , system of equations

p3a2:map("=",makelist(x,i,1,length(p2b)),p2b);

system:

p3a3:makelist(subst([p3a2[i]],p3a1),i,1,length(p3a2));

→seventh step :solve system "p3a3"

use :p2d2

p3a4:flatten(solve(p3a3,p2d2));

final result : numerical

→eight step :subst "p3a4" into "p3a1"

p4:subst(p3a4,p3a1);

→night step : error =exact-numerical

p4a:error = rhs(p1a2) -rhs(p4);

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