How to resolve the algorithm Numerical integration/Gauss-Legendre Quadrature step by step in the Axiom programming language

Published on 12 May 2024 09:40 PM
#Axiom

How to resolve the algorithm Numerical integration/Gauss-Legendre Quadrature step by step in the Axiom programming language

Table of Contents

Problem Statement

For this, we first need to calculate the nodes and the weights, but after we have them, we can reuse them for numerious integral evaluations, which greatly speeds up the calculation compared to more simple numerical integration methods.

Task description Similar to the task Numerical Integration, the task here is to calculate the definite integral of a function

f ( x )

{\displaystyle f(x)}

, but by applying an n-point Gauss-Legendre quadrature rule, as described here, for example. The input values should be an function f to integrate, the bounds of the integration interval a and b, and the number of gaussian evaluation points n. An reference implementation in Common Lisp is provided for comparison. To demonstrate the calculation, compute the weights and nodes for an 5-point quadrature rule and then use them to compute:

Let's start with the solution:

Step by Step solution about How to resolve the algorithm Numerical integration/Gauss-Legendre Quadrature step by step in the Axiom programming language

Source code in the axiom programming language

NNI ==> NonNegativeInteger
RECORD ==> Record(x : List Fraction Integer, w : List Fraction Integer)

gaussCoefficients(n : NNI, eps : Fraction Integer) : RECORD == 
  p := legendreP(n,z)
  q := n/2*D(p, z)*legendreP(subtractIfCan(n,1)::NNI, z)
  x := map(rhs,solve(p,eps))
  w := [subst(1/q, z=xi) for xi in x]
  [x,w]

gaussIntegrate(e : Expression Float, segbind : SegmentBinding(Float), n : NNI) : Float ==
  eps := 1/10^100
  u := gaussCoefficients(n,eps)
  interval := segment segbind
  var := variable segbind
  a := lo interval
  b := hi interval
  c := (a+b)/2
  h := (b-a)/2
  h*reduce(+,[wi*subst(e,var=c+xi*h) for xi in u.x for wi in u.w])

digits(50)
gaussIntegrate(4/(1+x^2), x=0..1, 20)

   (1)  3.1415926535_8979323846_2643379815_9534002592_872901276
                                                                  Type: Float
% - %pi

   (2)  - 0.3463549483_9378821092_475 E -26

You may also check:How to resolve the algorithm Gotchas step by step in the Raku programming language
You may also check:How to resolve the algorithm Continued fraction/Arithmetic/G(matrix ng, continued fraction n1, continued fraction n2) step by step in the Kotlin programming language
You may also check:How to resolve the algorithm Transliterate English text using the Greek alphabet step by step in the Raku programming language
You may also check:How to resolve the algorithm Read a file line by line step by step in the Vala programming language
You may also check:How to resolve the algorithm Power set step by step in the Déjà Vu programming language