feat: euler and picard-poly

euler-methods.py

@@ -0,0 +1,68 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.scale import FuncScale
+
+# f: (t, y) -> R, dy/dt = f
+def euler(f, a, b, y0, h):
+    xs = [a]
+    ys = [y0]
+    x = a
+    y = y0
+    while x < b:
+        y += h * f(x, y)
+        x += h
+        if abs(y) > threshold:
+            break
+        xs.append(x)
+        ys.append(y)
+    return xs, ys
+
+def improvedEuler(f, a, b, y0, h):
+    xs = [a]
+    ys = [y0]
+    
+    x = a
+    y = y0
+    while x <= b:
+        slope1 = f(x, y)
+        y_predictor = y + h * slope1
+        slope2 = f(x + h, y_predictor)
+        y = y + (h/2) * (slope1 + slope2)
+        x = x + h
+        if abs(y) > threshold:
+            break
+        xs.append(x)
+        ys.append(y)
+    return xs, ys
+
+f = lambda t, y: t**3 + y**3
+threshold = 100
+a = 0
+b = 2
+
+lims = [
+    (0, 1.4, 1.8, 1, threshold),
+    (1, 0.4, 0.6, 1, threshold),
+    (-1, 0.4, 0.6, -threshold, -1),
+]
+
+y0, xm, xM, ym, yM = lims[0]
+
+fig, ax= plt.subplots(1, 2)
+for h in [0.1, 0.01, 0.001, 0.0001]:
+    ax[0].plot(*euler(f, a, b, y0, h), label=f'$\\Delta t$={h}', linewidth=1)
+    T, Y = improvedEuler(f, a, b, y0, h)
+    ax[1].plot(T, Y, label=f'$\\Delta t$={h}', linewidth=1)
+    print(f"Max x: {T[-1]}")
+
+for axe in ax:
+    axe.legend()
+    axe.set_xlim(xm, xM)
+    axe.set_ylim(ym, yM)
+    #axe.set_yscale(FuncScale(axe, (lambda y: np.log(-y), lambda y: -np.exp(y))))
+
+ax[0].set_title("Euler Method")
+ax[1].set_title("Improved Euler Method")
+
+#fig.suptitle(r"Numerically solving $y'=y^3 + t^3$ starting at (0, 0)")
+plt.show()

picard-poly.py

@@ -0,0 +1,222 @@
+import numpy as np
+from numbers import Number
+import matplotlib.pyplot as plt
+
+class Polynomial:
+    """
+    Represents a polynomial and supports common mathematical operations.
+
+    The coefficients are stored in an array in order of increasing power,
+    e.g., `[c₀, c₁, c₂]` for the polynomial c₀ + c₁x + c₂x².
+    """
+
+    def __init__(self, coeff):
+        """Initializes the polynomial with a list, tuple, or numpy array of coefficients."""
+        # Trim trailing zeros that don't affect the polynomial's degree
+        self.coefficients = np.trim_zeros(np.asarray(coeff, dtype=float), 'b')
+        if self.coefficients.size == 0:
+            self.coefficients = np.array([0.0])
+
+
+    @property
+    def degree(self) -> int:
+        """Returns the degree of the polynomial."""
+        # The degree is the number of coefficients minus one.
+        # Returns -1 for the zero polynomial P(x) = 0.
+        if self.coefficients.size == 1 and self.coefficients[0] == 0:
+            return -1
+        return self.coefficients.size - 1
+
+    def evalf(self, x: float) -> float:
+        """
+        Evaluates the polynomial at a given value x using Horner's method.
+
+        For a polynomial P(x) = c₀ + c₁x + ... + cₙxⁿ, this is computed as
+        c₀ + x(c₁ + x(c₂ + ...)).
+        """
+        # The user's original implementation was correct.
+        rval = self.coefficients[-1]
+        for c in reversed(self.coefficients[:-1]):
+            rval = rval * x + c
+        return rval
+
+    def integrate(self) -> 'Polynomial':
+        """
+        Computes the indefinite integral of the polynomial, assuming the
+        constant of integration is zero.
+        """
+        # Corrected to use np.zeros for safe array initialization.
+        integrated_coeffs = np.zeros(self.coefficients.size + 1)
+        for i, c in enumerate(self.coefficients):
+            integrated_coeffs[i + 1] = c / (i + 1)
+        return Polynomial(integrated_coeffs)
+
+    def integratef(self, a: float, b: float) -> float:
+        """Computes the definite integral of the polynomial from a to b."""
+        # Corrected to pass the integration bounds to evalf.
+        integral = self.integrate()
+        return integral.evalf(b) - integral.evalf(a)
+
+    # --- Refinements and New Features ---
+
+    def __str__(self) -> str:
+        """Provides a user-friendly string representation of the polynomial."""
+        if self.degree == -1:
+            return "0"
+        
+        terms = []
+        # Iterate from the highest degree term to the lowest
+        for i, c in list(enumerate(self.coefficients)):
+            if np.isclose(c, 0) and i != 0:
+                continue
+
+            # Sign of the term
+            sign = ""
+            if len(terms) > 0:
+                sign = " + " if c > 0 else " - "
+            elif c < 0:
+                sign = "-"
+            
+            c = abs(c)
+
+            # Coefficient string
+            coeff_str = ""
+            if not np.isclose(c, 1) or i == 0:
+                coeff_str = f"{c:g}"
+
+            # Variable and power string
+            power_str = ""
+            if i > 0:
+                power_str = "x"
+                if i > 1:
+                    power_str += f"^{i}"
+            
+            # Combine coefficient and power
+            term = ""
+            if coeff_str and power_str:
+                term = coeff_str + "*" + power_str
+            else:
+                term = coeff_str or power_str
+
+            terms.append(sign + term)
+            
+        return "".join(terms)
+
+    def __repr__(self) -> str:
+        """Provides an unambiguous representation of the polynomial."""
+        return f"Polynomial({self.coefficients.tolist()})"
+
+    def __add__(self, other):
+        """Adds two polynomials or a polynomial and a scalar."""
+        if isinstance(other, Polynomial):
+            n1, n2 = self.coefficients.size, other.coefficients.size
+            if n1 > n2:
+                new_coeffs = self.coefficients.copy()
+                new_coeffs[:n2] += other.coefficients
+            else:
+                new_coeffs = other.coefficients.copy()
+                new_coeffs[:n1] += self.coefficients
+            return Polynomial(new_coeffs)
+        elif isinstance(other, Number):
+            new_coeffs = self.coefficients.copy()
+            new_coeffs[0] += other
+            return Polynomial(new_coeffs)
+        return NotImplemented
+
+    def __sub__(self, other):
+        """Subtracts two polynomials or a scalar from a polynomial."""
+        if isinstance(other, Polynomial):
+            n1, n2 = self.coefficients.size, other.coefficients.size
+            new_coeffs = np.zeros(max(n1, n2), dtype=float)
+            new_coeffs[:n1] = self.coefficients
+            new_coeffs[:n2] -= other.coefficients
+            return Polynomial(new_coeffs)
+        elif isinstance(other, Number):
+            new_coeffs = self.coefficients.copy()
+            new_coeffs[0] -= other
+            return Polynomial(new_coeffs)
+        return NotImplemented
+
+    def __mul__(self, other):
+        """Multiplies two polynomials or a polynomial and a scalar."""
+        if isinstance(other, Polynomial):
+            # Use convolution to multiply polynomial coefficients
+            return Polynomial(np.convolve(self.coefficients, other.coefficients))
+        elif isinstance(other, Number):
+            return Polynomial(self.coefficients * other)
+        return NotImplemented
+
+    def __pow__(self, n: int):
+        """Raises the polynomial to a non-negative integer power."""
+        if not isinstance(n, int) or n < 0:
+            raise ValueError("Exponent must be a non-negative integer.")
+        
+        # Use exponentiation by squaring for efficiency (O(log n) multiplications)
+        if n == 0:
+            return Polynomial([1]) # P^0 = 1
+        
+        result = Polynomial([1])
+        base = self
+        while n > 0:
+            if n % 2 == 1:
+                result *= base
+            base *= base
+            n //= 2
+        return result
+
+    def __radd__(self, other):
+        return self.__add__(other)
+
+    def __rsub__(self, other):
+        # other - self = -(self - other)
+        res = self.__sub__(other)
+        return res * -1
+
+    def __rmul__(self, other):
+        return self.__mul__(other)
+    
+    def __call__(self, t):
+        return self.evalf(t)
+    
+def picardIteration(t0, y0, num_iters = 20):
+    phi = [Polynomial([y0])]
+    for _ in range(num_iters): 
+        fty = Polynomial([0, 0, 0, 1]) + phi[-1]**3
+        intg = fty.integrate()
+        intg_result = intg - intg.evalf(t0)
+        phi.append(intg_result + phi[0])
+    return phi
+
+lims = [
+    (0, 2, 0, 100),
+    (0, 1, 1, 100),
+    (0, 1, -1, 100)
+]
+num_iter = 6
+
+t0, tM, y0, yM = lims[1]
+tm, ym = min(t0, tM), min(y0, yM)
+
+picard_result = picardIteration(t0, y0, num_iter)
+'''
+for idx, poly in enumerate(picard_result):
+    np_poly = np.vectorize(poly)
+    ts = np.linspace(t0, tM, 1000)
+    ys = np_poly(ts)
+    plt.plot(ts, ys)
+    print(idx, poly.degree, poly)
+'''
+poly = picard_result[-1]
+#np_poly = np.vectorize(poly)
+#ts = np.linspace(t0, tM, 1000)
+#ys = np_poly(ts)
+#plt.plot(ts, ys)
+#print(poly.degree, poly)
+print(poly.degree)
+plt.scatter(list(range(poly.degree + 1)), poly.coefficients, s=0.5)
+plt.yscale('log')
+#plt.xlim(t0, tM)
+#plt.ylim(ym, yM)
+plt.title(f"coefficients of $\\sum c_n t^n$ at y(0) = {y0} of iterations {num_iter}")
+plt.ylabel("$c_n$")
+plt.show()