Add UB4

UB2/QuadratureRule.cpp

@@ -5,13 +5,16 @@ QuadratureRule::QuadratureRule(const std::vector<double>& w, const std::vector<d
     : weights(w), points(x) {}
 
     // the integrate method takes a callable function f and integrates it over the domain R
-double QuadratureRule::integrate(const std::function<double(double)>& f) {
-    double integral = 0.0;
+double QuadratureRule::integrate(const std::function<double(double)>& f, double a, double b) {
+    double integral = 0;
+    double stepSize = (b-a)/2;
 
     // the integration is a weighted sum of function values at the specified points
     for (size_t i = 0; i < weights.size(); i++){
         integral += weights[i] * f(points[i]);
     }
     
+    integral *= stepSize;
+
     return integral;
 }
\ No newline at end of file

UB2/QuadratureRule.h

@@ -12,7 +12,7 @@ private:
 
 public:
     QuadratureRule(const std::vector<double>& w, const std::vector<double>& x);
-    double integrate(const std::function<double(double)>& f);
+    double integrate(const std::function<double(double)>& f, double a, double b);
 };
 
 #endif // QUADRATURERULE_H

UB2/main.cpp

@@ -4,12 +4,14 @@
 int main() {
     std::function<double(double)> func = [](double x) { return x * x; };
 
-    std::vector<double> weights = {1/3, 4/3, 1/3};
+    std::vector<double> weights = {1.0/3.0, 4.0/3.0, 1.0/3.0};
     std::vector<double> points = {-1,0,1};
 
     QuadratureRule quadrature(weights, points);
+    double a = -1;
+    double b = 1;
 
-    double result = quadrature.integrate(func);
+    double result = quadrature.integrate(func,a,b);
 
     std::cout << "Approximate integral of f(x) = x^2 over [-1, 1] is: " << result << std::endl;
 

UB2/polynomial_new.py

+import numpy as np
+
+class polynomial:
+    
+    def __init__(self, coeff):
+        self.coeff = coeff
+        self.deg = len(coeff) - 1
+        self._diff_coeff = np.zeros(self.deg + 1)
+        self._diff_coeff[:-1] = np.arange(1, self.deg + 1) * self.coeff[1:]
+        
+    def __call__(self, x):
+        if isinstance(x, float):
+            x = np.array([x])
+            X = self._evaluate_powers(x)
+            return (self.coeff @ X)[0]
+        else:
+            X = self._evaluate_powers(x)
+            return self.coeff @ X
+    
+    def diff_p(self, x):
+        X = self._evaluate_powers(x)
+        
+        return self._diff_coeff @ X
+        
+    def _evaluate_powers(self, x):
+        X = np.zeros((self.deg + 1, len(x)))
+        for ii in range(self.deg + 1):
+            X[ii, :] = x**ii
+        
+        return X
+
+class polynomial_new(polynomial):
+    
+    def __init__(self, coeff):
+        # This method calls the constructor of the parent class:
+        super(polynomial_new, self).__init__(coeff)
+        # Overwrite array of derivative coefficients, they are replaced by
+        # a matrix of coefficients for all (relevant) derivatives:
+        self._diff_coeff = self._eval_diff_coeff()
+        
+    
+    # Overwrite derivative method:
+    def diff_p(self, x, order=0):
+        
+        if order > self.deg:
+            return np.zeros(x.shape)
+        else:
+            # Evaluate all powers of x at all data sites:
+            X = self._evaluate_powers(x)
+
+            return self._diff_coeff[order, :] @ X
+    
+    # New function to generate derivative coefficients:
+    def _eval_diff_coeff(self):
+        # Compute all derivative coefficients:
+        diff_coeff = np.zeros((self.deg + 1, self.deg + 1))
+        diff_coeff[0, :] = self.coeff
+        for ii in range(1, self.deg + 1):
+            if ii == 1:
+                upper_ind = None
+            else:
+                upper_ind = -ii + 1
+            diff_coeff[ii, :-ii] = np.arange(1, self.deg + 2 - ii) * diff_coeff[ii-1, 1:upper_ind]
+        
+        return diff_coeff
\ No newline at end of file

UB2/quadrature_full.py

+import numpy as np
+
+class QuadratureRule:
+
+    def __init__(self, x, w):
+        self.x = x
+        self.w = w
+        self._n = self.x.shape[0]
+
+    def integrate(self, f):
+        # Evaluate f on all grid points:
+        fx = np.array([f(self.x[i]) for i in range(self._n)])
+
+        return np.dot(self.w, fx)
+    
+
+
+class LinearTransformed(QuadratureRule):
+
+    def integrate(self, f, T, dT):
+        # Evaluate f on transformed grid points:
+        fx = np.array([f(T(self.x[i])) for i in range(self._n)])
+        # Transform weights:
+        wT = np.array([self.w[i] * dT for i in range(self._n)])
+
+        return np.dot(wT, fx)
+
+
+class CompositeQuadrature:
+
+    def __init__(self, a, b, n, x, w):
+        """ Composite quadrature on interval [a, b] using n sub-
+         intervals, and a quadrature rule on reference interval [-1, 1].
+        
+        Parameters:
+        a, b:   lower and upper bounds of integration domain
+        n:      number of sub-intervals
+        x:      grid points for reference quadrature rule
+        w:      weights for reference quadrature rule
+        """
+        self.a = a
+        self.b = b
+        self.n = n
+        # Create sub-division of [a, b]:
+        self.h = (b - a) / n
+        self.y = np.array([a + i * self.h for i in range(self.n+1)])
+        # Mid-points of the sub-intervals:
+        self.c = 0.5 * (self.y[:-1] + self.y[1:])
+        # Create quadrature rule object:
+        self.QR = LinearTransformed(x, w)
+
+    def integrate(self, f):
+        I = 0.0
+        # Iterate over sub-intervals:
+        for i in range(self.n):
+            # Linear transformation:
+            T = lambda x: 0.5 * self.h * x + self.c[i]
+            I += self.QR.integrate(f, T, 0.5*self.h)
+        
+        return I
\ No newline at end of file

UB2/shell_script.sh

+g++ main.cpp QuadratureRule.cpp -o MyIntegrator
+./MyIntegrator

UB2/test_quad.py

+import numpy as np
+import quadrature
+import polynomial_new as poly
+
+""" Settings: """
+# Create an instance of the quadrature rule:
+x = np.array([-1.0, 0.0, 1.0])
+w = np.array([1.0/3, 4.0/3, 1.0/3])
+S = quadrature.QuadratureRule(x, w)
+
+# Boundary values of interval:
+a = -1.0
+b = 1.0
+# Number of test:
+ntest = 10
+
+""" Function to integrate polynomial: """
+def int_poly(P, a, b):
+    # Coefficients of the primitive function:
+    pc = P.coeff / np.arange(1, P.deg+2, dtype=float)
+    # Evaluate powers at both boundary points:
+    bx = np.array([b**i for i in range(1, P.deg+2)])
+    ax = np.array([a**i for i in range(1, P.deg+2)])
+
+    return (pc @ bx - pc @ ax)
+
+""" Main Test: """
+for ii in range(ntest):
+    # Generate random third order polynomial:
+    #c = np.array([1.0, 0.0, 0.0, 1.0])
+    c = np.random.randn(4)
+    print("Test %d: coeffs: "%ii, c)
+    P = poly.polynomial(c)
+    # Compute integral:
+    I = int_poly(P, a, b)
+    # Approximation with Simpson rule:
+    ISimp = S.integrate(P)
+    print("Test %d: diff = "%ii, np.abs(I - ISimp))
\ No newline at end of file

UB2/test_quad_composite.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+import quadrature_full
+
+""" Settings: """
+# Function to integrate:
+f = lambda x: np.arctan(x)
+# Anti-derivative:
+F = lambda x: x * np.arctan(x) - 0.5 * np.log(1 + x**2)
+
+# Interval:
+a = -1.0
+b = 3.0
+
+# Array of grid sizes:
+n_array = np.array([5, 10, 50, 100, 500, 1000, 5000, 10000])
+# Reference quadrature:
+x = np.array([-1.0, 0.0, 1.0])
+w = np.array([1.0/3, 4.0/3, 1.0/3])
+
+""" Demonstration: """
+# Compute reference:
+Iref = F(b) - F(a)
+# Convert grid sizes to step sizes:
+h_array = (b - a) / n_array
+# Apply composite quadrature:
+E = np.zeros(h_array.shape[0])
+for i in range(h_array.shape[0]):
+    n = n_array[i]
+    QC = quadrature_full.CompositeQuadrature(a, b, n, x, w)
+    I = QC.integrate(f)
+    E[i] = np.abs(I - Iref)
+
+""" Figure: """
+plt.figure(figsize=(6, 4))
+plt.plot(h_array, E, "o--", label="Error Simpson")
+plt.plot(h_array, h_array**4, "ro--", label=r"$h^4$")
+plt.xscale("log")
+plt.yscale("log")
+plt.xlabel(r"$h$")
+plt.ylim([1e-16, 1e-2])
+plt.legend(loc=2)
+plt.show()
+
+print(I, Iref)
\ No newline at end of file

UB2/test_quad_transformed.py

+import numpy as np
+import quadrature
+import polynomial_new as poly
+
+""" Settings: """
+# Parameters of the quadrature rule on reference interval:
+x = np.array([-1.0, 0.0, 1.0])
+w = np.array([1.0/3, 4.0/3, 1.0/3])
+
+# Define linear transformation of [-1, 1] to [1, 5]:
+a = 1.0
+b = 5.0
+T = lambda x: 2 * x + 3
+dT = 2.0
+# Number of tests:
+ntest = 10
+
+""" Function to integrate polynomial: """
+def int_poly(P, a, b):
+    # Coefficients of the primitive function:
+    pc = P.coeff / np.arange(1, P.deg+2, dtype=float)
+    # Evaluate powers at both boundary points:
+    bx = np.array([b**i for i in range(1, P.deg+2)])
+    ax = np.array([a**i for i in range(1, P.deg+2)])
+
+    return (pc @ bx - pc @ ax)
+
+
+""" Main Test: """
+S = quadrature.LinearTransformed(x, w)
+
+for ii in range(ntest):
+    # Generate random third order polynomial:
+    #c = np.array([1.0, 0.0, 0.0, 1.0])
+    c = np.random.randn(4)
+    print("Test %d: coeffs: "%ii, c)
+    P = poly.polynomial(c)
+    # Compute integral:
+    I = int_poly(P, a, b)
+    # Approximation with Simpson rule:
+    ISimp = S.integrate(P, T, dT)
+    print("Test %d: diff = "%ii, np.abs(I - ISimp))
\ No newline at end of file

UB3/random_walk.py

+import numpy as np
+import random
+import matplotlib.pyplot as plt
+
+from run_rw2d import _update_figure
+
+class RandomWalk2D:
+    def __init__(self, dx: float, N: int, x0: np.ndarray, dt: float):
+        """
+        :param x0: vector of initial positions
+        :param N: population size
+        """
+        assert x0.shape == (N,2), "x0 must have the shape (N,2)"
+        self.positions = x0
+        self.time = [0]
+        self.dt = dt
+        self.N = N
+        self.dx = dx
+    
+    def simulate(self, K: int = 1):
+        """
+        param K: number of discrete steps
+        """
+
+        choices = [(-1, 0), (1, 0), (0, 1), (0, -1)]
+
+        for _ in range(K):
+            self.time.append(self.time[-1]+ self.dt)
+            directions = np.array([random.choice(choices) for _ in range(self.N)])* self.dx
+            self.positions += directions
+
+
+### testing ###
+N=100
+K=20
+dt=4e-2
+dx=4e-1   
+RW = RandomWalk2D(dx=dx, N=N, x0 = np.zeros((N,2)), dt=dt)
+
+### visualisation ###
+
+# Axes limits for figure:
+xmin = -8.0
+xmax = 8.0
+X = RW.positions
+fig = plt.figure(figsize=(6, 4))
+fig = _update_figure(fig, X, 0.0, xmin, xmax)
+for _ in range(1, K+1):
+    # Update positions:
+    RW.simulate()
+    X = RW.positions
+    print(f'{X=}')
+    t = RW.time[-1]
+    # Update figure:
+    fig = _update_figure(fig, X, t, xmin, xmax)
+    plt.pause(0.25)
+
+plt.show()

UB3/run_rw2d.py

+import numpy as np
+import matplotlib.pyplot as plt
+
+# import random_walk_2d as rw
+
+""" Settings: """
+# Mesh sizes:
+dx = 0.4
+dt = 0.04
+# Population size:
+N = 1000
+# Initial conditions:
+x0 = np.zeros((2, N))
+t0 = 0.0
+# Number of discrete time steps:
+K = 20
+# Axes limits for figure:
+xmin = -8.0
+xmax = 8.0
+
+""" Function to update the figure: """
+def _update_figure(fig, X, t, xmin, xmax):
+    fig.clear()
+    plt.scatter(X[:, 0], X[:, 1])
+    plt.title("Population at time t = %.3f"%t)
+    plt.xlim([xmin, xmax])
+    plt.ylim([xmin, xmax])
+
+    return fig
+
+# """ Instantiate the class: """
+# RW = rw.RandomWalk2D(dx, dt, N, x0)
+
+# """ Simulation: """
+# X = RW.positions()
+# fig = plt.figure(figsize=(6, 4))
+# fig = _update_figure(fig, X, 0.0, xmin, xmax)
+
+# for ii in range(1, K+1):
+#     # Update positions:
+#     RW.simulate()
+#     X = RW.positions()
+#     t = RW.time()
+#     # Update figure:
+#     fig = _update_figure(fig, X, t, xmin, xmax)
+#     plt.pause(0.25)
+
+# plt.show()

UB4/UB4.py

+import numpy as np
+import matplotlib.pyplot  as plt
+
+# Ex1: Numerical integrators 
+
+def implicit_euler(f, y0: np.ndarray, t: float, dt: float) -> np.ndarray:
+    """
+    A-stable
+    :param f: function to be integrated
+    :param y0: initial value
+    :param t: time interval
+    :param dt: time step
+    """
+    no_of_steps = int(t // dt)
+    y = np.zeros((no_of_steps, *y0.shape))
+    y[0] = y0
+    
+    for i in range(1, no_of_steps):
+        y[i] = y[i-1] + f(y[i-1]) * dt
+        y[i] = y[i-1] + f(y[i]) * dt
+
+    return y
+
+def explicit_euler(f, y0: np.ndarray, t: float, dt: float) -> np.ndarray:
+    """
+    :param f: function to be integrated
+    :param y0: initial value
+    :param t: time interval
+    :param dt: time step
+    """
+    no_of_steps = int(t // dt)
+    y = np.zeros((no_of_steps, *y0.shape))
+    y[0] = y0
+
+    for i in range(1, no_of_steps):
+        y[i] = y[i-1] + f(y[i-1]) * dt
+    
+    return y
+
+def implicit_midpoint(f, y0: np.ndarray, t: float, dt: float) -> np.ndarray:
+    """
+    Not L-stable, doesn't decay properly - oscillation
+    :param f: function to be integrated
+    :param y0: initial value
+    :param t: time interval
+    :param dt: time step
+    """
+    no_of_steps = int(t // dt)
+    y = np.zeros((no_of_steps, *y0.shape))
+    y[0] = y0
+    
+    for i in range(1, no_of_steps):
+        y[i] = y[i-1] + dt * f(y[i-1] + dt/2 * f(y[i-1]))
+        y[i] = y[i-1] + dt * f(1/2 *  (y[i-1] + y[i]))
+    
+    return y
+
+# Ex2: dy/dt = lambda * y, y(0) = 1, lambda = -1
+y0 = np.array([1])
+lamda = -1
+t = 100
+dt = 1e-2
+f = lambda y: lamda * y
+y = [implicit_euler(f, y0, t, dt), explicit_euler(f, y0, t, dt), implicit_midpoint(f, y0, t ,dt)]
+integrators = ['implicit euler', 'explicit euler', 'implicit midpoint']
+
+# for i in range(3):
+#     plt.figure()
+#     plt.plot(y[i])
+#     plt.title(integrators[i])
+#     plt.show()
+
+# Ex3: 
+def f(y: np.ndarray, k: float = 1) -> np.ndarray:
+    assert y.shape[0] == 3, 'y has the wrong dimension. It should be 3'
+
+    A = np.array([[-1, 1, 0],
+                    [1, -1-k, k],
+                    [0, k, -k]])
+    
+    return np.dot(A, y)
+
+dt = [1, 1e-1, 1e-2, 1e-3]
+t = [10,100]
+y0 = np.array([1,0,0])
+
+# short time
+fig = plt.figure()
+for dt_i in dt:
+    y = implicit_euler(f, y0, t[0], dt_i)
+    plt.plot(y[:,0], '--',label = 'State 1')
+    plt.plot(y[:,1], 'x', label = 'State 2')
+    plt.plot(y[:,2], '>', label = 'State 3')
+    plt.legend()
+    plt.title(f'{integrators[0]}, dt = {dt_i}')
+    plt.show()
+
+
+
+
+
+
+fig = plt.figure()
+for dt_i in dt:
+    y = explicit_euler(f, y0, t[0], dt_i)
+    plt.plot(y[:,0], '--', label = 'State 1')
+    plt.plot(y[:,1], 'x', label = 'State 2')
+    plt.plot(y[:,2], '>', label = 'State 3')
+    plt.legend()
+    plt.title(f'{integrators[1]}, dt = {dt_i}')
+    plt.show()
+
+fig = plt.figure()
+for dt_i in dt:
+    y = implicit_midpoint(f, y0, t[0], dt_i)
+    plt.plot(y[:,0], '--', label = 'State 1')
+    plt.plot(y[:,1], 'x', label = 'State 2')
+    plt.plot(y[:,2], '>', label = 'State 3')
+    plt.legend()
+    plt.title(f'{integrators[2]}, dt = {dt_i}')
+    plt.show()
+
+

shell_script.sh

-g++ main.cpp QuadratureRule.cpp -o MyApp
-./MyIntegrator