diff --git a/UB2/MyIntegrator b/UB2/MyIntegrator
index 026b9ee91cf3c24edc6d2e0cd1b732abd64b7a49..50c1a570180b960467c9c23fed2122ee6a925bda 100755
Binary files a/UB2/MyIntegrator and b/UB2/MyIntegrator differ
diff --git a/UB2/QuadratureRule.cpp b/UB2/QuadratureRule.cpp
index 3063ce3d9a36a75a91b0b606ab3125b1dc9e508a..7591a8de2bd0424d5a297974c0128388f7b49198 100644
--- a/UB2/QuadratureRule.cpp
+++ b/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
diff --git a/UB2/QuadratureRule.h b/UB2/QuadratureRule.h
index 134ae1a31e428c3e3ea5ad0d2d0cc6895af68258..f90d8b3aee369b142399ad1f2ee524b9196214a7 100644
--- a/UB2/QuadratureRule.h
+++ b/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
diff --git a/UB2/main.cpp b/UB2/main.cpp
index 9709cf5797eced7e31b51c738abd8758663e3415..1501401880bf4647a9779f9f1f768f8305fdd267 100644
--- a/UB2/main.cpp
+++ b/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;
diff --git a/UB2/polynomial_new.py b/UB2/polynomial_new.py
new file mode 100644
index 0000000000000000000000000000000000000000..310787899e587834d1161308dc92e696dd662f36
--- /dev/null
+++ b/UB2/polynomial_new.py
@@ -0,0 +1,65 @@
+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
diff --git a/UB2/quadrature_full.py b/UB2/quadrature_full.py
new file mode 100644
index 0000000000000000000000000000000000000000..dbe04a0cc3c15fd54cf31f270551f393f9d0b117
--- /dev/null
+++ b/UB2/quadrature_full.py
@@ -0,0 +1,60 @@
+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
diff --git a/UB2/shell_script.sh b/UB2/shell_script.sh
new file mode 100644
index 0000000000000000000000000000000000000000..1c3085330c44f04ad39e51260175ce65e4d72ce1
--- /dev/null
+++ b/UB2/shell_script.sh
@@ -0,0 +1,2 @@
+g++ main.cpp QuadratureRule.cpp -o MyIntegrator
+./MyIntegrator
diff --git a/UB2/test_quad.py b/UB2/test_quad.py
new file mode 100644
index 0000000000000000000000000000000000000000..280f47f37c4af467d9151397b3a8b96fc56d38cf
--- /dev/null
+++ b/UB2/test_quad.py
@@ -0,0 +1,38 @@
+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
diff --git a/UB2/test_quad_composite.py b/UB2/test_quad_composite.py
new file mode 100644
index 0000000000000000000000000000000000000000..e56026ee37c19772f4a6baa44ec56cdfb8f11b1f
--- /dev/null
+++ b/UB2/test_quad_composite.py
@@ -0,0 +1,46 @@
+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
diff --git a/UB2/test_quad_transformed.py b/UB2/test_quad_transformed.py
new file mode 100644
index 0000000000000000000000000000000000000000..c03aadb91203c5915f5b5594b969432813a01b27
--- /dev/null
+++ b/UB2/test_quad_transformed.py
@@ -0,0 +1,42 @@
+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
diff --git a/UB3/__pycache__/run_rw2d.cpython-310.pyc b/UB3/__pycache__/run_rw2d.cpython-310.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..eb602d2a62b12a374abbab7324c6114c71e11aed
Binary files /dev/null and b/UB3/__pycache__/run_rw2d.cpython-310.pyc differ
diff --git a/UB3/random_walk.py b/UB3/random_walk.py
new file mode 100644
index 0000000000000000000000000000000000000000..3fd6fbb4963e74929f84de07bde458548508441b
--- /dev/null
+++ b/UB3/random_walk.py
@@ -0,0 +1,58 @@
+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()
diff --git a/UB3/run_rw2d.py b/UB3/run_rw2d.py
new file mode 100644
index 0000000000000000000000000000000000000000..24ac93b7b5bea6b0187e3f9f3e3aef2fcdd8d441
--- /dev/null
+++ b/UB3/run_rw2d.py
@@ -0,0 +1,48 @@
+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()
diff --git a/UB4/UB4.py b/UB4/UB4.py
new file mode 100644
index 0000000000000000000000000000000000000000..402e74d170d879f7d2f022f3ee6d13c364ce385a
--- /dev/null
+++ b/UB4/UB4.py
@@ -0,0 +1,123 @@
+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()
+
+
diff --git a/shell_script.sh b/shell_script.sh
deleted file mode 100644
index 6d1b172c5a5264eadcd945b2b999bc3b9984156d..0000000000000000000000000000000000000000
--- a/shell_script.sh
+++ /dev/null
@@ -1,2 +0,0 @@
-g++ main.cpp QuadratureRule.cpp -o MyApp
-./MyIntegrator