Skip to content

Instantly share code, notes, and snippets.

@HASPIMA
Last active November 3, 2018 04:27
Show Gist options
  • Select an option

  • Save HASPIMA/a77a469cdca52b89f95ba1b8d5f8af09 to your computer and use it in GitHub Desktop.

Select an option

Save HASPIMA/a77a469cdca52b89f95ba1b8d5f8af09 to your computer and use it in GitHub Desktop.
This software was made for my integral calculus class. Its main purpose is to use the Gauss–Legendre quadrature for make an approximation to integral.
# /usr/bin/python3
from math import e, sqrt
# coefficients has the respective values for each n
coefficients = (((0.0, 2.0), ()), # ((xi, wi))
((1 / sqrt(3), 1), (-1 / sqrt(3), 1)),
((0, 8 / 9), (sqrt(3 / 5), 5 / 9), (-sqrt(3 / 5), 5 / 9)),
(
(sqrt((3 - 2 * sqrt(6 / 5)) / 7), (18 + sqrt(30)) / 36),
(-sqrt((3 - 2 * sqrt(6 / 5)) / 7), (18 + sqrt(30)) / 36),
(sqrt((3 + 2 * sqrt(6 / 5)) / 7), (18 - sqrt(30)) / 36),
(-sqrt((3 + 2 * sqrt(6 / 5)) / 7), (18 - sqrt(30)) / 36)
),
(
(0, 128 / 225),
(sqrt(5 - 2 * sqrt(10 / 7)) / 3, (322 + 13 * sqrt(70)) / 900),
(-sqrt(5 - 2 * sqrt(10 / 7)) / 3, (322 + 13 * sqrt(70)) / 900),
(sqrt(5 + 2 * sqrt(10 / 7)) / 3, (322 - 13 * sqrt(70)) / 900),
(-sqrt(5 + 2 * sqrt(10 / 7)) / 3, (322 - 13 * sqrt(70)) / 900)
)
)
# function of testing
def f(x: float) -> float:
return e ** -x ** 2
def quadrature(function_, a: float, b: float, nodes: int = 3):
if nodes < 1 or nodes > 5:
return None
c1, c2 = (b - a) / 2, (a + b) / 2
amt = 0
for i in range(nodes):
amt += coefficients[nodes - 1][i][1] * function_(c1 * coefficients[nodes - 1][i][0] + c2)
return c1 * amt
def main():
from sys import stdin, stdout
a, b, n = map(float, stdin.readline().rstrip().split())
n = int(n)
stdout.write(str(quadrature(f, a, b, n)) + '\n')
if __name__ == '__main__':
main()
@HASPIMA
Copy link
Author

HASPIMA commented Nov 3, 2018

Now it works!

You only have to input a, b and nodes respectively. Remember to change the parameter of the line 47 for your own function that you want to integrate with the Gauss–Legendre quadrature.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment