From 7a2a843285e86b4822ebc53cef2013fd5879edce Mon Sep 17 00:00:00 2001 From: "Tammo.Weber" Date: Fri, 6 Feb 2026 14:10:48 +0100 Subject: [PATCH] kleinere Optimierungen --- GHA_triaxial/gha2_num.py | 331 +++++++++++++++------------------------ runge_kutta.py | 82 +++++++++- 2 files changed, 203 insertions(+), 210 deletions(-) diff --git a/GHA_triaxial/gha2_num.py b/GHA_triaxial/gha2_num.py index b912ba1..a700fe0 100644 --- a/GHA_triaxial/gha2_num.py +++ b/GHA_triaxial/gha2_num.py @@ -1,24 +1,19 @@ import numpy as np from ellipsoide import EllipsoidTriaxial -import runge_kutta as rk +from runge_kutta import rk4, rk4_step, rk4_end, rk4_integral import GHA_triaxial.numeric_examples_karney as ne_karney import GHA_triaxial.numeric_examples_panou as ne_panou import winkelumrechnungen as wu from typing import Tuple from numpy.typing import NDArray import ausgaben as aus -from utils_angle import cot, wrap_to_pi +from utils_angle import cot, arccot, wrap_to_pi -def arccot(x): - x = np.asarray(x) - a = np.arctan2(1.0, x) - return np.where(x < 0.0, a - np.pi, a) - -def normalize_alpha_0_pi(alpha): - if alpha < 0.0: - alpha += np.pi - return alpha +def norm_a(a): + if a < 0.0: + a += np.pi + return a def sph_azimuth(beta1, lam1, beta2, lam2): dlam = wrap_to_pi(lam2 - lam1) @@ -32,17 +27,29 @@ def sph_azimuth(beta1, lam1, beta2, lam2): # Panou 2013 def gha2_num( ell: EllipsoidTriaxial, + beta_0: float, + lamb_0: float, beta_1: float, lamb_1: float, - beta_2: float, - lamb_2: float, n: int = 16000, - epsilon: float = 10 ** -12, + epsilon: float = 10**-12, iter_max: int = 30, all_points: bool = False, ) -> Tuple[float, float, float] | Tuple[float, float, float, NDArray, NDArray]: + """ + :param ell: Ellipsoid + :param beta_0: Beta Punkt 0 + :param lamb_0: Lambda Punkt 0 + :param beta_1: Beta Punkt 1 + :param lamb_1: Lambda Punkt 1 + :param n: Anzahl Schritte + :param epsilon: Genauigkeit + :param iter_max: Maximale Iterationen + :param all_points: Ausgabe aller Punkte + :return: Azimut Startpunkt, Azumit Zielpunkt, Strecke + """ - + # Berechnung Koeffizienten, Gaußschen Fundamentalgrößen 1. Ordnung sowie deren Ableitungen def BETA_LAMBDA(beta, lamb): BETA = (ell.ay ** 2 * np.sin(beta) ** 2 + ell.b ** 2 * np.cos(beta) ** 2) / ( ell.Ex ** 2 - ell.Ey ** 2 * np.sin(beta) ** 2 @@ -123,6 +130,7 @@ def gha2_num( G_lamb_lamb, ) + # Berechnung der ODE Koeffizienten für Fall 1 (lambda_0 != lambda_1) def p_coef(beta, lamb): ( BETA, @@ -161,6 +169,7 @@ def gha2_num( return (BETA, LAMBDA, E, G, p_3, p_2, p_1, p_0, p_33, p_22, p_11, p_00) + # Berechnung der ODE Koeffizienten für Fall 2 (lambda_0 == lambda_1) def q_coef(beta, lamb): ( BETA, @@ -197,69 +206,7 @@ def gha2_num( ) q_00 = 0.5 * ((E_lamb_lamb * G - E_lamb * G_lamb) / (G**2)) - return (BETA, LAMBDA, E, G, q_3, q_2, q_1, q_0, q_33, q_22, q_11, q_00) - - def rk4_last(f, t0, y0, dt, N): - h = dt / N - t = t0 - y = np.array(y0, dtype=float, copy=True) - for _ in range(N): - k1 = f(t, y) - k2 = f(t + 0.5 * h, y + 0.5 * h * k1) - k3 = f(t + 0.5 * h, y + 0.5 * h * k2) - k4 = f(t + h, y + h * k3) - y = y + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4) - t = t + h - return t, y - - def rk4_last_with_integral(f, t0, y0, dt, N, integrand_at): - - h = dt / N - habs = abs(h) - t = t0 - y = np.array(y0, dtype=float, copy=True) - - if N % 2 == 0: - # Simpson streaming - f0 = integrand_at(t, y) - odd_sum = 0.0 - even_sum = 0.0 - - for i in range(1, N + 1): - k1 = f(t, y) - k2 = f(t + 0.5 * h, y + 0.5 * h * k1) - k3 = f(t + 0.5 * h, y + 0.5 * h * k2) - k4 = f(t + h, y + h * k3) - y = y + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4) - t = t + h - - fi = integrand_at(t, y) - if i == N: - fN = fi - elif i % 2 == 1: - odd_sum += fi - else: - even_sum += fi - - S = f0 + fN + 4.0 * odd_sum + 2.0 * even_sum - s = (habs / 3.0) * S - return t, y, s - - # Trapez streaming - f_prev = integrand_at(t, y) - acc = 0.0 - for _ in range(N): - k1 = f(t, y) - k2 = f(t + 0.5 * h, y + 0.5 * h * k1) - k3 = f(t + 0.5 * h, y + 0.5 * h * k2) - k4 = f(t + h, y + h * k3) - y = y + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4) - t = t + h - f_cur = integrand_at(t, y) - acc += 0.5 * (f_prev + f_cur) - f_prev = f_cur - s = habs * acc - return t, y, s + return BETA, LAMBDA, E, G, q_3, q_2, q_1, q_0, q_33, q_22, q_11, q_00 def integrand_lambda(lamb, y): beta = y[0] @@ -273,42 +220,41 @@ def gha2_num( (_, _, E, G, *_) = BETA_LAMBDA(beta, lamb) return np.sqrt(E + G * lamb_p**2) - if lamb_1 != lamb_2: - N = n - dlamb = lamb_2 - lamb_1 + # Fall 1 (lambda_0 != lambda_1) + if abs(lamb_1 - lamb_0) >= 1e-15: + N = int(n) + dlamb = float(lamb_1 - lamb_0) - def buildODElamb(): - def ODE(lamb, v): - beta, beta_p, X3, X4 = v - (_, _, _, _, p_3, p_2, p_1, p_0, p_33, p_22, p_11, p_00) = p_coef(beta, lamb) + beta0 = float(beta_0) + lamb0 = float(lamb_0) + beta1 = float(beta_1) + lamb1 = float(lamb_1) - dbeta = beta_p - dbeta_p = p_3 * beta_p**3 + p_2 * beta_p**2 + p_1 * beta_p + p_0 - dX3 = X4 - dX4 = (p_33 * beta_p**3 + p_22 * beta_p**2 + p_11 * beta_p + p_00) * X3 + ( - 3 * p_3 * beta_p**2 + 2 * p_2 * beta_p + p_1 - ) * X4 - return np.array([dbeta, dbeta_p, dX3, dX4], dtype=float) + def ode_lamb(lamb, v): + beta, beta_p, X3, X4 = v + (_, _, _, _, p_3, p_2, p_1, p_0, p_33, p_22, p_11, p_00) = p_coef(beta, lamb) - return ODE - - ode_lamb = buildODElamb() - - alpha0_sph = sph_azimuth(beta_1, lamb_1, beta_2, lamb_2) - (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1) - beta_p0_sph = np.sqrt(G1 / E1) * cot(alpha0_sph) if abs(dlamb) >= 1e-15 else 0.0 + dbeta = beta_p + dbeta_p = p_3 * beta_p**3 + p_2 * beta_p**2 + p_1 * beta_p + p_0 + dX3 = X4 + dX4 = (p_33 * beta_p**3 + p_22 * beta_p**2 + p_11 * beta_p + p_00) * X3 + (3*p_3*beta_p**2 + 2*p_2*beta_p + p_1) * X4 + return np.array([dbeta, dbeta_p, dX3, dX4], dtype=float) + alpha0_sph = sph_azimuth(beta0, lamb0, beta1, lamb1) + (_, _, E0, G0, *_) = BETA_LAMBDA(beta0, lamb0) + beta_p0_sph = np.sqrt(G0 / E0) * cot(alpha0_sph) N_newton = min(N, 4000) def solve_newton(beta_p0_init: float): beta_p0 = float(beta_p0_init) for _ in range(iter_max): - startwerte = np.array([beta_1, beta_p0, 0.0, 1.0], dtype=float) - _, y_end = rk4_last(ode_lamb, lamb_1, startwerte, dlamb, N_newton) - beta_end, _, X3_end, _ = y_end + v0 = np.array([beta0, beta_p0, 0.0, 1.0], dtype=float) + _, y_end = rk4_end(ode_lamb, lamb0, v0, dlamb, N_newton) + + beta_end, _, X3_end, _ = y_end + delta = beta_end - beta1 - delta = beta_end - beta_2 if abs(delta) < epsilon: return True, beta_p0 @@ -316,194 +262,161 @@ def gha2_num( return False, None step = delta / X3_end - max_step = 0.5 - if abs(step) > max_step: - step = np.sign(step) * max_step - beta_p0 = beta_p0 - step + step = np.clip(step, -0.5, 0.5) + beta_p0 -= step return False, None ok, beta_p0_sol = solve_newton(beta_p0_sph) if not ok: - candidates = [-beta_p0_sph, 0.5 * beta_p0_sph, 2.0 * beta_p0_sph] N_quick = min(N, 2000) best = None - for g in candidates: - ok_g, beta_p0_sol_g = solve_newton(g) + ok_g, sol = solve_newton(g) if not ok_g: continue - - startwerte_g = np.array([beta_1, beta_p0_sol_g, 0.0, 1.0], dtype=float) - _, _, s_quick = rk4_last_with_integral( - ode_lamb, lamb_1, startwerte_g, dlamb, N_quick, integrand_lambda - ) - + v0_g = np.array([beta0, sol, 0.0, 1.0], dtype=float) + _, _, s_quick = rk4_integral(ode_lamb, lamb0, v0_g, dlamb, N_quick, integrand_lambda) if (best is None) or (s_quick < best[0]): - best = (s_quick, beta_p0_sol_g) - + best = (s_quick, sol) if best is None: - raise RuntimeError("Keine Startwert-Variante konvergiert.") - + raise RuntimeError("Keine Startwert-Variante konvergiert (lambda-Fall).") beta_p0_sol = best[1] - beta_0 = beta_p0_sol - startwerte_final = np.array([beta_1, beta_0, 0.0, 1.0], dtype=float) + beta_p0 = float(beta_p0_sol) + v0_final = np.array([beta0, beta_p0, 0.0, 1.0], dtype=float) if all_points: - lamb_list, states = rk.rk4(ode_lamb, lamb_1, startwerte_final, dlamb, N, False) - + lamb_list, states = rk4(ode_lamb, lamb0, v0_final, dlamb, N, False) lamb_arr = np.array(lamb_list, dtype=float) beta_arr = np.array([st[0] for st in states], dtype=float) beta_p_arr = np.array([st[1] for st in states], dtype=float) - (_, _, E1, G1, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) - (_, _, E2, G2, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1]) + (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) + (_, _, E_end, G_end, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1]) - alpha_1 = arccot(np.sqrt(E1 / G1) * beta_p_arr[0]) - alpha_2 = arccot(np.sqrt(E2 / G2) * beta_p_arr[-1]) + alpha_1 = norm_a(arccot(np.sqrt(E_start / G_start) * beta_p_arr[0])) + alpha_2 = norm_a(arccot(np.sqrt(E_end / G_end) * beta_p_arr[-1])) - alpha_1 = normalize_alpha_0_pi(float(alpha_1)) - alpha_2 = normalize_alpha_0_pi(float(alpha_2)) - - # Distanz s aus Arrays (Simpson/Trapz) + # Distanz aus Arrays integrand = np.zeros(N + 1, dtype=float) for i in range(N + 1): (_, _, Ei, Gi, *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i]) - integrand[i] = np.sqrt(Ei * beta_p_arr[i] ** 2 + Gi) + integrand[i] = np.sqrt(Ei * beta_p_arr[i]**2 + Gi) h = abs(dlamb) / N if N % 2 == 0: - S = integrand[0] + integrand[-1] + 4.0 * np.sum(integrand[1:-1:2]) + 2.0 * np.sum(integrand[2:-1:2]) - s = h / 3.0 * S + S = integrand[0] + integrand[-1] + 4.0*np.sum(integrand[1:-1:2]) + 2.0*np.sum(integrand[2:-1:2]) + s = h/3.0 * S else: s = np.trapz(integrand, dx=h) - return alpha_1, alpha_2, s, beta_arr, lamb_arr + return float(alpha_1), float(alpha_2), float(s), beta_arr, lamb_arr - _, y_end, s = rk4_last_with_integral(ode_lamb, lamb_1, startwerte_final, dlamb, N, integrand_lambda) + _, y_end, s = rk4_integral(ode_lamb, lamb0, v0_final, dlamb, N, integrand_lambda) beta_end, beta_p_end, _, _ = y_end - (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1) - (_, _, E2, G2, *_) = BETA_LAMBDA(beta_2, lamb_2) + (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta0, lamb0) + (_, _, E_end, G_end, *_) = BETA_LAMBDA(beta1, lamb1) - alpha_1 = arccot(np.sqrt(E1 / G1) * beta_0) - alpha_2 = arccot(np.sqrt(E2 / G2) * beta_p_end) + alpha_1 = norm_a(arccot(np.sqrt(E_start / G_start) * beta_p0)) + alpha_2 = norm_a(arccot(np.sqrt(E_end / G_end) * beta_p_end)) - alpha_1 = normalize_alpha_0_pi(float(alpha_1)) - alpha_2 = normalize_alpha_0_pi(float(alpha_2)) + return float(alpha_1), float(alpha_2), float(s) - return alpha_1, alpha_2, s - - N = n - dbeta = beta_2 - beta_1 + # Fall 2 (lambda_0 == lambda_1) + N = int(n) + dbeta = float(beta_1 - beta_0) if abs(dbeta) < 1e-15: if all_points: return 0.0, 0.0, 0.0, np.array([]), np.array([]) return 0.0, 0.0, 0.0 - # ODE-System (lambda, lambda', Y3, Y4) in Abhängigkeit von beta - def buildODEbeta(): - def ODE(beta, v): - lamb, lamb_p, Y3, Y4 = v - (_, _, _, _, q_3, q_2, q_1, q_0, q_33, q_22, q_11, q_00) = q_coef(beta, lamb) + beta0 = float(beta_0) + lamb0 = float(lamb_0) + beta1 = float(beta_1) + lamb1 = float(lamb_1) - dlamb = lamb_p - dlamb_p = q_3 * lamb_p**3 + q_2 * lamb_p**2 + q_1 * lamb_p + q_0 - dY3 = Y4 - dY4 = (q_33 * lamb_p**3 + q_22 * lamb_p**2 + q_11 * lamb_p + q_00) * Y3 + ( - 3 * q_3 * lamb_p**2 + 2 * q_2 * lamb_p + q_1 - ) * Y4 + def ode_beta(beta, v): + lamb, lamb_p, Y3, Y4 = v + (_, _, _, _, q_3, q_2, q_1, q_0, q_33, q_22, q_11, q_00) = q_coef(beta, lamb) - return np.array([dlamb, dlamb_p, dY3, dY4], dtype=float) + dlamb = lamb_p + dlamb_p = q_3*lamb_p**3 + q_2*lamb_p**2 + q_1*lamb_p + q_0 + dY3 = Y4 + dY4 = (q_33*lamb_p**3 + q_22*lamb_p**2 + q_11*lamb_p + q_00)*Y3 + (3*q_3*lamb_p**2 + 2*q_2*lamb_p + q_1)*Y4 + return np.array([dlamb, dlamb_p, dY3, dY4], dtype=float) - return ODE - - ode_beta = buildODEbeta() - - # Newton auf lambda'_0 - lamb_0 = 0.0 + lamb_p0 = 0.0 for _ in range(iter_max): - startwerte = np.array([lamb_1, lamb_0, 0.0, 1.0], dtype=float) - beta_list, states = rk.rk4(ode_beta, beta_1, startwerte, dbeta, N, False) + v0 = np.array([lamb0, lamb_p0, 0.0, 1.0], dtype=float) + _, y_end = rk4_end(ode_beta, beta0, v0, dbeta, N) - lamb_end, lamb_p_end, Y3_end, _ = states[-1] - delta = lamb_end - lamb_2 + lamb_end, _, Y3_end, _ = y_end + delta = lamb_end - lamb1 if abs(delta) < epsilon: break if abs(Y3_end) < 1e-20: - raise RuntimeError("Abbruch (Ableitung ~ 0).") + raise RuntimeError("Abbruch (Ableitung ~ 0) im beta-Fall.") step = delta / Y3_end - max_step = 1.0 - if abs(step) > max_step: - step = np.sign(step) * max_step - lamb_0 = lamb_0 - step + step = np.clip(step, -1.0, 1.0) + lamb_p0 -= step - startwerte_final = np.array([lamb_1, lamb_0, 0.0, 1.0], dtype=float) + v0_final = np.array([lamb0, lamb_p0, 0.0, 1.0], dtype=float) if all_points: - beta_list, states = rk.rk4(ode_beta, beta_1, startwerte_final, dbeta, N, False) - + beta_list, states = rk4(ode_beta, beta0, v0_final, dbeta, N, False) beta_arr = np.array(beta_list, dtype=float) lamb_arr = np.array([st[0] for st in states], dtype=float) lamb_p_arr = np.array([st[1] for st in states], dtype=float) - # Azimute - (BETA1, LAMBDA1, _, _, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) - (BETA2, LAMBDA2, _, _, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1]) + (BETA_s, LAMBDA_s, _, _, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) + (BETA_e, LAMBDA_e, _, _, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1]) - alpha_1 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA1 / BETA1) * lamb_p_arr[0]) - alpha_2 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA2 / BETA2) * lamb_p_arr[-1]) + alpha_1 = norm_a((np.pi/2.0) - arccot(np.sqrt(LAMBDA_s / BETA_s) * lamb_p_arr[0])) + alpha_2 = norm_a((np.pi/2.0) - arccot(np.sqrt(LAMBDA_e / BETA_e) * lamb_p_arr[-1])) - # optionaler Quadrantenfix (robust) - alpha_1 = normalize_alpha_0_pi(float(alpha_1)) - alpha_2 = normalize_alpha_0_pi(float(alpha_2)) - - # Distanz integrand = np.zeros(N + 1, dtype=float) for i in range(N + 1): (_, _, Ei, Gi, *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i]) - integrand[i] = np.sqrt(Ei + Gi * lamb_p_arr[i] ** 2) + integrand[i] = np.sqrt(Ei + Gi * lamb_p_arr[i]**2) h = abs(dbeta) / N if N % 2 == 0: - S = integrand[0] + integrand[-1] + 4.0 * np.sum(integrand[1:-1:2]) + 2.0 * np.sum(integrand[2:-1:2]) - s = h / 3.0 * S + S = integrand[0] + integrand[-1] + 4.0*np.sum(integrand[1:-1:2]) + 2.0*np.sum(integrand[2:-1:2]) + s = h/3.0 * S else: s = np.trapz(integrand, dx=h) - return alpha_1, alpha_2, s, beta_arr, lamb_arr + return float(alpha_1), float(alpha_2), float(s), beta_arr, lamb_arr - # all_points == False: streaming Integral - _, y_end, s = rk4_last_with_integral(ode_beta, beta_1, startwerte_final, dbeta, N, integrand_beta) + _, y_end, s = rk4_integral(ode_beta, beta0, v0_final, dbeta, N, integrand_beta) lamb_end, lamb_p_end, _, _ = y_end - (BETA1, LAMBDA1, _, _, *_) = BETA_LAMBDA(beta_1, lamb_1) - (BETA2, LAMBDA2, _, _, *_) = BETA_LAMBDA(beta_2, lamb_2) + (BETA_s, LAMBDA_s, _, _, *_) = BETA_LAMBDA(beta0, lamb0) + (BETA_e, LAMBDA_e, _, _, *_) = BETA_LAMBDA(beta1, lamb1) - alpha_1 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA1 / BETA1) * lamb_0) - alpha_2 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA2 / BETA2) * lamb_p_end) - - alpha_1 = normalize_alpha_0_pi(float(alpha_1)) - alpha_2 = normalize_alpha_0_pi(float(alpha_2)) - - return alpha_1, alpha_2, s + alpha_1 = norm_a((np.pi/2.0) - arccot(np.sqrt(LAMBDA_s / BETA_s) * lamb_p0)) + alpha_2 = norm_a((np.pi/2.0) - arccot(np.sqrt(LAMBDA_e / BETA_e) * lamb_p_end)) + return float(alpha_1), float(alpha_2), float(s) if __name__ == "__main__": - ell = EllipsoidTriaxial.init_name("BursaSima1980round") - beta1 = np.deg2rad(75) - lamb1 = np.deg2rad(-90) - beta2 = np.deg2rad(75) - lamb2 = np.deg2rad(66) - a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=5000) - print(aus.gms("a1", a1, 4)) + # ell = EllipsoidTriaxial.init_name("BursaSima1980round") + # beta1 = np.deg2rad(75) + # lamb1 = np.deg2rad(-90) + # beta2 = np.deg2rad(75) + # lamb2 = np.deg2rad(66) + # a0, a1, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=5000) + # print(aus.gms("a0", a0, 4)) + # print(aus.gms("a1", a1, 4)) + # print("s: ", s) # # print(aus.gms("a2", a2, 4)) # # print(s) # cart1 = ell.para2cart(0, 0) diff --git a/runge_kutta.py b/runge_kutta.py index 4ab408e..e786507 100644 --- a/runge_kutta.py +++ b/runge_kutta.py @@ -40,4 +40,84 @@ def rk4_step(ode, t: float, v: np.ndarray, h: float) -> np.ndarray: k2 = ode(t + 0.5 * h, v + 0.5 * h * k1) k3 = ode(t + 0.5 * h, v + 0.5 * h * k2) k4 = ode(t + h, v + h * k3) - return v + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4) \ No newline at end of file + return v + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4) + +def rk4_end(ode, t0: float, v0: np.ndarray, weite: float, schritte: int, fein: bool = False): + h = weite / schritte + t = float(t0) + v = np.array(v0, dtype=float, copy=True) + + for _ in range(schritte): + if not fein: + v_next = rk4_step(ode, t, v, h) + else: + v_grob = rk4_step(ode, t, v, h) + v_half = rk4_step(ode, t, v, 0.5 * h) + v_fein = rk4_step(ode, t + 0.5 * h, v_half, 0.5 * h) + v_next = v_fein + (v_fein - v_grob) / 15.0 + + t += h + v = v_next + + return t, v + +# RK4 mit Simpson bzw. Trapez +def rk4_integral( ode, t0: float, v0: np.ndarray, weite: float, schritte: int, integrand_at, fein: bool = False, simpson: bool = True, ): + + h = weite / schritte + habs = abs(h) + + t = float(t0) + v = np.array(v0, dtype=float, copy=True) + + if simpson and (schritte % 2 == 0): + f0 = float(integrand_at(t, v)) + odd_sum = 0.0 + even_sum = 0.0 + fN = None + + for i in range(1, schritte + 1): + if not fein: + v_next = rk4_step(ode, t, v, h) + else: + v_grob = rk4_step(ode, t, v, h) + v_half = rk4_step(ode, t, v, 0.5 * h) + v_fein = rk4_step(ode, t + 0.5 * h, v_half, 0.5 * h) + v_next = v_fein + (v_fein - v_grob) / 15.0 + + t += h + v = v_next + + fi = float(integrand_at(t, v)) + if i == schritte: + fN = fi + elif i % 2 == 1: + odd_sum += fi + else: + even_sum += fi + + S = f0 + fN + 4.0 * odd_sum + 2.0 * even_sum + s = (habs / 3.0) * S + return t, v, s + + f_prev = float(integrand_at(t, v)) + acc = 0.0 + + for _ in range(schritte): + if not fein: + v_next = rk4_step(ode, t, v, h) + else: + v_grob = rk4_step(ode, t, v, h) + v_half = rk4_step(ode, t, v, 0.5 * h) + v_fein = rk4_step(ode, t + 0.5 * h, v_half, 0.5 * h) + v_next = v_fein + (v_fein - v_grob) / 15.0 + + t += h + v = v_next + + f_cur = float(integrand_at(t, v)) + acc += 0.5 * (f_prev + f_cur) + f_prev = f_cur + + s = habs * acc + return t, v, s \ No newline at end of file