<sxh python> #! /usr/bin/env python # -*- coding: utf-8 -*- # vim:fenc=utf-8 # # Copyright © 2017 putanowr <putanowr@foo> # # Distributed under terms of the MIT license.
“”“ Newton-Cotes numerical integration formulas ”“”
class NewtonCotes:
def __init__(self, degree): self.degree = degree def nodes(self, a, b): n = self.degree+1 dx = (b-a)/n return [a+i*dx for i in range(n)]
def weights(self, a, b): return ((b-a)*x for x in self.ref_weights)
class TrapezoidQuadrature(NewtonCotes):
ref_weights = (0.5, 0.5)
def __init__(self):
super().__init__(1)
class SimpsonQuadrature(NewtonCotes):
ref_weights = (x/6.0 for x in (1, 4, 1))
def __init__(self):
super().__init__(2)
if name == 'main':
qr = TrapezoidQuadrature()
a = -3;
b = 5;
if abs(sum(qr.weights(a,b)) - (b-a)) < 1.e-8:
print("Tesing weights: OK")
else:
print("Testing weights: FAILED")
print("TrapezoidQuadrature degree: ", qr.degree)
qr = SimpsonQuadrature()
if abs(sum(qr.weights(a,b)) - (b-a)) < 1.e-8:
print("Tesing weights: OK")
else:
print("Testing weights: FAILED")
print("SimpsonQuadrature degree: ", qr.degree)
</sxh>
<sxh python> #! /usr/bin/env python # -*- coding: utf-8 -*- # vim:fenc=utf-8 # # Copyright © 2017 putanowr <putanowr@foo> # # Distributed under terms of the MIT license.
“”“ Calculate integral of scalar function in 1D ”“” import sys import os import itertools import newtoncotes
class Integrand:
"""Represents integrand as 1D scalar function """ def __init__(self, expression): self.expression = expression
def evaluate(self, x): return eval(self.expression)
def __call__(self, x): return self.evaluate(x)
class Mesh:
"""Represents on dimensional mesh """ def __init__(self): self.nodes = list()
def load(self, filename):
"""Read from file a list of nodes"""
with open(filename, 'r') as f:
for l in f:
self.nodes.append(tuple(float(x) for x in l.split()))
def nelem(self): return len(self.nodes)-1
class MeshIntegrator:
"""Calculate integral over a domain discretised by a mesh """
def integrate(self, mesh, integrand):
"""Integrate function fun using numerical integration on given mesh
"""
quadrature = newtoncotes.TrapezoidQuadrature();
integral = 0.0
for i in range(mesh.nelem()):
integral += self.integrate_element(i, mesh, integrand, quadrature)
return integral
def make_mesh_fun(self, mesh, fun):
xcoord = [ node[0] for node in mesh.nodes ] discrete_fun = [ fun(x) for x in xcoord ] return discrete_fun
def integrate_element(self, i, mesh, integrand, quadrature):
x1 = mesh.nodes[i][0]
x2 = mesh.nodes[i+1][0]
qr = quadrature
integral = 0.0
for (x, w) in itertools.zip_longest(qr.nodes(x1, x2), qr.weights(x1,x2)):
integral += w * integrand(x);
return integral
def parse_command_line(argv):
"""Parse command line"""
argc = len(sys.argv)
if len(sys.argv) < 2:
print("Mesh file name must be given")
sys.exit(22)
meshfile = sys.argv[1]
if len(sys.argv) > 2:
fun = sys.argv[2]
else:
print("Using default function f(x) = x^2")
fun = "x**2"
return [meshfile, fun]
def main():
[meshfile, fun] = parse_command_line(sys.argv)
mesh = Mesh()
mesh.load(meshfile)
integrand = Integrand(fun)
integrator = MeshIntegrator()
integral = integrator.integrate(mesh, integrand)
print("Integral of %s is %g" %(fun, integral))
if name == 'main':
main()
</sxh>