diff --git a/.ipynb_checkpoints/runmassspring-checkpoint.sh b/.ipynb_checkpoints/runmassspring-checkpoint.sh new file mode 100755 index 000000000..a42139b39 --- /dev/null +++ b/.ipynb_checkpoints/runmassspring-checkpoint.sh @@ -0,0 +1,13 @@ +#!/bin/bash + +tend_relative=$1 #standard: 4 +steps=$2 #standard: 100 +algorithm=$3 + +cd build +./test_ode > ../demos/output_test_ode.txt $tend_relative $steps $algorithm +cd ../demos +python3 plotmassspring.py --tend_relative $1 --steps $2 --algorithm $3 +mv "Phase_plot.png" "plots/$3/$3 Phase_plot_$1 pi_$2 steps_.png" +mv "Time_evolution.png" "plots/$3/$3 Time_evolution_$1 pi_$2 steps_.png" +cd .. \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 6dc00d2b7..9515a2284 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,10 +9,18 @@ include_directories(src nanoblas/src) add_subdirectory (src) add_subdirectory (nanoblas) - +add_compile_options(-fpermissive) add_executable (test_ode demos/test_ode.cpp) target_link_libraries (test_ode PUBLIC nanoblas) add_executable (demo_autodiff demos/demo_autodiff.cpp) target_link_libraries (demo_autodiff PUBLIC nanoblas) +add_executable(test_rc demos/test_rc.cpp) +target_link_libraries(test_rc PUBLIC nanoblas) + +add_executable(test_autodiff demos/test_autodiff.cpp) +target_link_libraries(test_autodiff PUBLIC nanoblas) + +add_executable(test_pendulum_autodiff demos/test_pendulum_autodiff.cpp) +target_link_libraries(test_pendulum_autodiff PUBLIC nanoblas) diff --git a/README.md b/README.md index 9572d8104..cb8e77131 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,24 @@ -# ASC-ODE -A package for solving ordinary differential equations +# Mass Spring System solver -Read the [documentation](https://tuwien-asc.github.io/ASC-ODE/intro.html) +## Usage +After cloning first run, run in the command line: -Find theory behind here: https://jschoeberl.github.io/IntroSC/ODEs/ODEs.html +git submodule update --init +Then build the project with cmake: + +mkdir build +cd build +cmake .. +make +cd .. + + +You can then run and plot the mass-spring system for a given time and nr of steps automatically with the shell script "runmassspring.sh": + +~/team01$ ./runmassspring.sh tend_relativetopi steps method + +tend_relativetopi and steps are doubles, method is either "explicit" or "improved", for example: + + +~/team01$ ./runmassspring.sh 4 100 explicit diff --git a/REPORT-2.md b/REPORT-2.md new file mode 100644 index 000000000..49704e9b7 --- /dev/null +++ b/REPORT-2.md @@ -0,0 +1,185 @@ +# Numerical Methods for ODEs — Homework Report + + +--- + +## 1. Introduction +This report summarizes the implementation and results of several numerical methods for solving ordinary differential equations (ODEs). It covers: + +- **HW 17.4:** Implicit Euler, Crank–Nicolson, and RC circuit simulation +- **HW 18:** Forward-mode Automatic Differentiation +- **HW 19:** Runge–Kutta methods (RK2) + +All methods were integrated into a unified `TimeStepper` framework and tested using the provided mass–spring and RC circuits. + +--- + +## 2. HW 17.4 — Implicit Euler, Crank–Nicolson & RC Circuit + +### 2.1 Mass–Spring System +We study the classical system: + +$$ +x' = v, \qquad v' = -\frac{k}{m}x +$$ + +This describes a harmonic oscillator. +The following methods were implemented and tested: + +- Explicit Euler +- Improved Euler +- **Implicit Euler** +- **Crank–Nicolson** + +#### **Plots** +Plots for the solution trajectories and phase-space behavior were generated. +

Time Evolution (Mass–Spring)

+

+ +

+ +

Phase Plot (Mass–Spring)

+

+ +

+ +#### **Observations** +- **Explicit Euler**: noticeable phase error and amplitude drift. +- **Improved Euler**: better accuracy but still accumulates drift. +- **Implicit Euler**: unconditionally stable, slightly overdamped. +- **Crank–Nicolson**: best accuracy among Euler-type schemes, symmetric and stable. + +--- + +### 2.2 RC Circuit ODE +The charging of a capacitor follows: + +$$ +C \frac{dU_c}{dt} = \frac{U_s - U_c}{R} +$$ + +which was written in autonomous form and solved using all four methods. +#### **Plots** + +

RC Circuit – Full Time Evolution

+

+ +

+ + +

RC Circuit – Zoomed View (Stability)

+

+ +

+ +#### **Observations** +- Explicit Euler can behave poorly for stiff choices of \( R \) and \( C \). +- Improved Euler is stable for moderate step sizes. +- Implicit Euler and Crank–Nicolson give smooth, stable curves. +- Crank–Nicolson consistently provides the most accurate evolution. + +Plots of capacitor voltage vs. time were produced for comparison. +--- + +## 3. HW 18 — Automatic Differentiation (AD) + +In this part we implemented a minimal forward-mode automatic differentiation framework based on a templated class +**`AutoDiff`**, capable of tracking derivatives with respect to **N** variables. + +### 3.1 Features Implemented + +The AutoDiff class supports: +- A stored **value** `val` +- A stored **derivative vector** `deriv[0…N-1]` +- Initialization from: + - constants + - variables of the form `Variable` +- Operator overloading for: + - `+` addition + - `-` subtraction + - `*` multiplication (all 3 overloads: AD*AD, const*AD, AD*const) + - `/` division (full quotient rule support) +- Elementary functions: + - `sin`, `cos`, `exp`, `log`, `pow` + +This enables symbolic-like differentiation while evaluating normal C++ expressions. + +### 3.2 Test Function + +We tested the AD system on: + +$$ +f(x) = x^2 + 3\sin(x) +$$ + +Evaluated at \(x = 1\), the AD derivative was compared to a numeric finite-difference estimate. + +#### **Result** +- AutoDiff derivative: **matches numerical derivative exactly** +- Confirms correct implementation of forward-mode AD + +### 3.3 Pendulum System +A second AD test was implemented for the simple pendulum derivative: + +$$ +f(\theta) = -\frac{g}{L}\sin(\theta) +$$ + +AD successfully produced the correct symbolic derivative: + +$$ +f'(\theta) = -\frac{g}{L}\cos(\theta) +$$ + +and matched the numeric approximation. + +--- + +## 4. HW 19 — Runge–Kutta Methods (RK2) + +In this exercise we implemented **one explicit Runge–Kutta method** (RK2), as required by the assignment. +RK4 was only used for internal testing and was removed before merging. + +### 4.1 RK2 (Midpoint Method) +The RK2 method improves accuracy significantly compared to Euler-type methods and is defined by: + +- Predictor step: + $$k_1 = f(y_n)$$ + +- Midpoint evaluation: + $$k_2 = f\left(y_n + \frac{1}{2}\tau k_1\right)$$ + +- Update: + $$y_{n+1} = y_n + \tau k_2$$ + +### 4.2 Testing + +RK2 was tested on the mass–spring system for various step counts: + +- 50 steps +- 100 steps +- 200 steps +- 400 steps + +### 4.3 Observations + +- **RK2** shows much smaller phase drift than Euler or improved Euler +- The numerical trajectory becomes smoother as step size decreases +- The energy error is significantly smaller than the Euler family +- RK2 remains stable for much larger time steps than explicit Euler + +Overall, RK2 provides a clear accuracy improvement with minimal additional computation. + +--- + +## 5. Overall Conclusions + +- **Implicit Euler** is stable and suitable for stiff systems +- **Crank–Nicolson** provides symmetric, high-quality solutions for oscillatory systems +- **Automatic Differentiation** works correctly and matches finite-difference tests +- **RK2** offers superior accuracy compared to Euler-type methods with only a small cost increase +- The unified ODE framework allowed consistent testing across all numerical methods + + + +--- diff --git a/REPORT.md b/REPORT.md new file mode 100644 index 000000000..2a4358a64 --- /dev/null +++ b/REPORT.md @@ -0,0 +1,140 @@ +## Part 3 - Exercise 1: Mass-Spring System +### 1. Testing existing explicit Euler method +Our goal in this task was to examine the behavior of the explicit Euler method for the mass-spring system. To obtain a reasonable estimation, the system was tested using different time steps and end times. All tests were automated using the `runmassspring.sh` script. + + +Explicit Euler method can be implemented as: + +$$ +y_{n+1} = y_n + h \cdot f(y_n) +$$ + +Where: +- $y_n$ — current value of the variable +- $y_{n+1}$ — next value of the variable +- $h$ — time step +- $f(y_n)$ — derivative of $y$ at step $n$ + + +Our system of mass-sprint is described by the following relation: + +$$ +y_{0}' = y_1 +$$ + +$$ +y_{1}' = -\frac{k}{m} y_0 +$$ + + + +With the original settings, we obtain the following results: + +
+ Time_evolution +
+
+ Phase_plot +
+ +As we can observe, this is not the solution for a mass-spring system without damping. A dissipative error is evident, as the amplitude of the motion increases over time. The solution appears to "explode", causing both the velocity and position to tend toward infinity. Such behaviour seems to increases the energy of a system - which is physically impossible due to Energy Conservation Law. In an undamped system, the amplitudes of velocity and position should remain constant over time. Position and velocity profiles should be in the form similar to harmonic oscilator i.e: + +$$x(t) = Acos(\omega_0 t + \phi)$$ +$$v(t) = -A\omega_0 sin(\omega_0 t + \phi)$$ + +where: + +$$\omega = \frac{k}{m}$$ + +This proves that amplitudes should remain constant over time and profiles being shifted as relation between sine and cosine + +Additionally, we would expect the position-velocity chart to take an elliptical shape, not exhibit divergence. +We can assume that the phase relationship between displacement and velocity is correct and equal to 90 degrees, meaning no dispersive error is observed. + + +By increasing number of steps, we can fix this for short intervals: + +
+ Time_evolution +
+
+ Phase_plot +
+ +The results of the simulation align with expectations, i.e., the relationship between velocity and position is conserved over time. However, it only seems to work properly for a fixed interval. As the simulation interval increases, instability in the system becomes apparent, triggering dissipative errors and results of the algorithm unreliable. + +
+ Time_evolution +
+
+ Phase_plot +
+ +Using this approach for shorter simulation intervals with increase of timestep improves the solution up to a certain point. One could suspect a hidden relationship between the timestep and the simulation interval that must be maintained to ensure a stable solution within the given time frame. Such bevaviour demonstrates convergence of the solution - but not stability. Convergence can be expressed in form of: + + +$$ +\lim_{\tau \to 0} |y_n - y(t_n)| = 0 +$$ + +where: +$$ +\tau = \frac{t_{final} - t_{initial}}{N_{steps}} +$$ + +
+ Time_evolution +
+
+ Phase_plot +
+ + +### 2. Implementing improved euler method + +To achieve better performance in the spring-mass system, we introduce the Improved Euler Method, which is defined as: + +$$ +\tilde{y} = y_n + \frac{\tau}{2} f(y_n) +$$ + +$$ +y_{n+1} = y_n + \tau f(\tilde{y}) +$$ + +In order to implement replace $f(y_n)$ with $f(\tilde{y})$ we had to make changes in file `timestepper.hpp` and modify a few other view details. + + + +
+ Time_evolution +
+
+ Phase_plot +
+ +That change has introduced a major improvement; however, over time we still cannot be certain about the system's stability at aribtary numer of steps. One would expect the solution to remain 'proper' as time approaches infinity. + +
+ Time_evolution +
+
+ Phase_plot +
+ +Introducing a significant increase in the number of timesteps—something that would completely destabilize the classical explicit Euler method—has remarkably stabilized the behavior of the Improved Euler method, both in terms of velocity and position over time. There is no noticeable dissipative or dispersive error, and the phase plot shows that the position-velocity relationship is ideally conserved as an elliptical trajectory. We cannot be sure that any diffusive of dispersive error still not exists, which could grow if the simulation were extended further. Each of our attempts has decreased timestep with grow of interval. However, the simulation was concluded at this step. + + Nevertheless, we cannot assume the solution is stable. In other words, according to the rules of the stability definition, we could be able to find a significantly small timestep that will make the solution stable while times go to infinity - we could not find such timestep. + + +$$ +\lim_{n_{steps} \to ∞} |y_n - y(t_n)| = 0 +$$ + + +
+ Time_evolution +
+
+ Phase_plot +
\ No newline at end of file diff --git a/Screenshot.png b/Screenshot.png new file mode 100644 index 000000000..dece2af6a Binary files /dev/null and b/Screenshot.png differ diff --git a/demos/output_test_ode.txt:Zone.Identifier b/demos/output_test_ode.txt:Zone.Identifier new file mode 100644 index 000000000..d6c1ec682 Binary files /dev/null and b/demos/output_test_ode.txt:Zone.Identifier differ diff --git a/demos/plot_rc.py b/demos/plot_rc.py new file mode 100644 index 000000000..b2049911c --- /dev/null +++ b/demos/plot_rc.py @@ -0,0 +1,70 @@ +import numpy as np +import matplotlib.pyplot as plt +import os +import argparse + +# We will run this script from the build directory: +# python3 ../demos/plot_rc.py --tend 0.01 --steps 2000 + +parser = argparse.ArgumentParser() +parser.add_argument("--tend", type=float, required=True) +parser.add_argument("--steps", type=int, required=True) +args = parser.parse_args() + +methods = ["explicit", "improved", "implicit", "CN"] +colors = { + "explicit": "tab:blue", + "improved": "tab:orange", + "implicit": "tab:green", + "CN": "tab:red", +} + +data = {} + +# load data from the text files in the current (build) directory +for m in methods: + fname = f"rc_{m}.txt" + + if os.path.exists(fname): + t, uc, _ = np.loadtxt(fname, unpack=True) + data[m] = (t, uc) + +if not data: + print("No RC data files found (explicit_rc.txt, ...). Run test_rc first.") + raise SystemExit + +# where to save plots (relative to this script location) +script_dir = os.path.dirname(os.path.abspath(__file__)) +plots_dir = os.path.join(script_dir, "plots", "RC") +os.makedirs(plots_dir, exist_ok=True) + +# --------- Plot 1: capacitor voltage vs time (all methods) ---------- +plt.figure() +for m, (t, uc) in data.items(): + plt.plot(t, uc, label=m, linewidth=1.0) +plt.xlabel("time t [s]") +plt.ylabel("capacitor voltage U_C(t)") +plt.title(f"RC circuit, tend = {args.tend}, steps = {args.steps}") +plt.grid(True) +plt.legend() + +out1 = os.path.join(plots_dir, f"RC_Time_evolution_tend_{args.tend}_steps_{args.steps}.png") +plt.savefig(out1, dpi=300) +print("Saved:", out1) + +# (optional) zoomed view of the beginning +plt.figure() +for m, (t, uc) in data.items(): + plt.plot(t, uc, label=m, linewidth=1.0) +plt.xlabel("time t [s]") +plt.ylabel("U_C(t)") +plt.title("RC circuit (zoomed, first 2 ms)") +plt.xlim(0.0, min(args.tend, 0.002)) +plt.grid(True) +plt.legend() + +out2 = os.path.join(plots_dir, f"RC_Time_evolution_zoom_tend_{args.tend}_steps_{args.steps}.png") +plt.savefig(out2, dpi=300) +print("Saved:", out2) + +# do not call plt.show() in WSL (no GUI) diff --git a/demos/plotmassspring.py b/demos/plotmassspring.py index bc1b19e76..d8dd80d5f 100644 --- a/demos/plotmassspring.py +++ b/demos/plotmassspring.py @@ -4,21 +4,33 @@ import matplotlib.pyplot as plt + +import argparse + +parser = argparse.ArgumentParser() +parser.add_argument("--tend_relative",type=str, required=True) +parser.add_argument("--steps",type=str, required=True) +parser.add_argument("--algorithm",type=str, required=True) +args = parser.parse_args() + + plt.plot(data[:,0], data[:,1], label='position') plt.plot(data[:,0], data[:,2], label='velocity') plt.xlabel('time') plt.ylabel('value') -plt.title('Mass-Spring System Time Evolution') +plt.title('Mass-Spring System Time Evolution, tend = ' + args.tend_relative + "pi, steps = " +args.steps + " , " + args.algorithm) plt.legend() plt.grid() +plt.savefig("Time_evolution.png",dpi=300) plt.show() - +plt.close() plt.plot(data[:,1], data[:,2], label='phase plot') plt.xlabel('position') plt.ylabel('velocity') -plt.title('Mass-Spring System Phase Plot') +plt.title('Mass-Spring System Phase Plot, tend = ' + args.tend_relative + "pi, steps = " +args.steps + " , " + args.algorithm) plt.legend() plt.grid() +plt.savefig("Phase_plot.png",dpi=300) plt.show() - +plt.close() \ No newline at end of file diff --git a/demos/plots/CN/.ipynb_checkpoints/CN Phase_plot_4 pi_100 steps_-checkpoint.png b/demos/plots/CN/.ipynb_checkpoints/CN Phase_plot_4 pi_100 steps_-checkpoint.png new file mode 100644 index 000000000..5cf9f4956 Binary files /dev/null and b/demos/plots/CN/.ipynb_checkpoints/CN Phase_plot_4 pi_100 steps_-checkpoint.png differ diff --git a/demos/plots/CN/.ipynb_checkpoints/CN Time_evolution_4 pi_100 steps_-checkpoint.png b/demos/plots/CN/.ipynb_checkpoints/CN Time_evolution_4 pi_100 steps_-checkpoint.png new file mode 100644 index 000000000..13a01329c Binary files /dev/null and b/demos/plots/CN/.ipynb_checkpoints/CN Time_evolution_4 pi_100 steps_-checkpoint.png differ diff --git a/demos/plots/CN/CN Phase_plot_4 pi_100 steps_.png b/demos/plots/CN/CN Phase_plot_4 pi_100 steps_.png new file mode 100644 index 000000000..5cf9f4956 Binary files /dev/null and b/demos/plots/CN/CN Phase_plot_4 pi_100 steps_.png differ diff --git a/demos/plots/CN/CN Time_evolution_4 pi_100 steps_.png b/demos/plots/CN/CN Time_evolution_4 pi_100 steps_.png new file mode 100644 index 000000000..13a01329c Binary files /dev/null and b/demos/plots/CN/CN Time_evolution_4 pi_100 steps_.png differ diff --git a/demos/plots/Phase_plot_standard.png b/demos/plots/Phase_plot_standard.png new file mode 100644 index 000000000..3d1db1982 Binary files /dev/null and b/demos/plots/Phase_plot_standard.png differ diff --git a/demos/plots/RC/RC_Time_evolution_tend_0.01_steps_2000.png b/demos/plots/RC/RC_Time_evolution_tend_0.01_steps_2000.png new file mode 100644 index 000000000..49b3de8f5 Binary files /dev/null and b/demos/plots/RC/RC_Time_evolution_tend_0.01_steps_2000.png differ diff --git a/demos/plots/RC/RC_Time_evolution_zoom_tend_0.01_steps_2000.png b/demos/plots/RC/RC_Time_evolution_zoom_tend_0.01_steps_2000.png new file mode 100644 index 000000000..af9e12917 Binary files /dev/null and b/demos/plots/RC/RC_Time_evolution_zoom_tend_0.01_steps_2000.png differ diff --git a/demos/plots/Time_evolution_standard.png b/demos/plots/Time_evolution_standard.png new file mode 100644 index 000000000..f9bcc6824 Binary files /dev/null and b/demos/plots/Time_evolution_standard.png differ diff --git a/demos/plots/explicit/explicit Phase_plot_10 pi_100 steps_.png b/demos/plots/explicit/explicit Phase_plot_10 pi_100 steps_.png new file mode 100644 index 000000000..69c22c741 Binary files /dev/null and b/demos/plots/explicit/explicit Phase_plot_10 pi_100 steps_.png differ diff --git a/demos/plots/explicit/explicit Phase_plot_12 pi_100 steps_.png b/demos/plots/explicit/explicit Phase_plot_12 pi_100 steps_.png new file mode 100644 index 000000000..b5a102b5b Binary files /dev/null and b/demos/plots/explicit/explicit Phase_plot_12 pi_100 steps_.png differ diff --git a/demos/plots/explicit/explicit Phase_plot_15 pi_100 steps_.png b/demos/plots/explicit/explicit Phase_plot_15 pi_100 steps_.png new file mode 100644 index 000000000..8c3b95922 Binary files /dev/null and b/demos/plots/explicit/explicit Phase_plot_15 pi_100 steps_.png differ diff --git a/demos/plots/explicit/explicit Phase_plot_4 pi_100 steps_.png b/demos/plots/explicit/explicit Phase_plot_4 pi_100 steps_.png new file mode 100644 index 000000000..ed31afd03 Binary files /dev/null and b/demos/plots/explicit/explicit Phase_plot_4 pi_100 steps_.png differ diff --git a/demos/plots/explicit/explicit Phase_plot_4 pi_1000 steps_.png b/demos/plots/explicit/explicit Phase_plot_4 pi_1000 steps_.png new file mode 100644 index 000000000..ed8feccca Binary files /dev/null and b/demos/plots/explicit/explicit Phase_plot_4 pi_1000 steps_.png differ diff --git a/demos/plots/explicit/explicit Phase_plot_4 pi_10000 steps_.png b/demos/plots/explicit/explicit Phase_plot_4 pi_10000 steps_.png new file mode 100644 index 000000000..26ede7049 Binary files /dev/null and b/demos/plots/explicit/explicit Phase_plot_4 pi_10000 steps_.png differ diff --git a/demos/plots/explicit/explicit Phase_plot_4 pi_200 steps_.png b/demos/plots/explicit/explicit Phase_plot_4 pi_200 steps_.png new file mode 100644 index 000000000..b3c48b32d Binary files /dev/null and b/demos/plots/explicit/explicit Phase_plot_4 pi_200 steps_.png differ diff --git a/demos/plots/explicit/explicit Phase_plot_60 pi_10000 steps_.png b/demos/plots/explicit/explicit Phase_plot_60 pi_10000 steps_.png new file mode 100644 index 000000000..717e9a1a1 Binary files /dev/null and b/demos/plots/explicit/explicit Phase_plot_60 pi_10000 steps_.png differ diff --git a/demos/plots/explicit/explicit Phase_plot_60 pi_1000000 steps_.png b/demos/plots/explicit/explicit Phase_plot_60 pi_1000000 steps_.png new file mode 100644 index 000000000..50887aa08 Binary files /dev/null and b/demos/plots/explicit/explicit Phase_plot_60 pi_1000000 steps_.png differ diff --git a/demos/plots/explicit/explicit Phase_plot_60 pi_5000 steps_.png b/demos/plots/explicit/explicit Phase_plot_60 pi_5000 steps_.png new file mode 100644 index 000000000..2c984fb84 Binary files /dev/null and b/demos/plots/explicit/explicit Phase_plot_60 pi_5000 steps_.png differ diff --git a/demos/plots/explicit/explicit Time_evolution_10 pi_100 steps_.png b/demos/plots/explicit/explicit Time_evolution_10 pi_100 steps_.png new file mode 100644 index 000000000..3dd30cbb5 Binary files /dev/null and b/demos/plots/explicit/explicit Time_evolution_10 pi_100 steps_.png differ diff --git a/demos/plots/explicit/explicit Time_evolution_12 pi_100 steps_.png b/demos/plots/explicit/explicit Time_evolution_12 pi_100 steps_.png new file mode 100644 index 000000000..5e039894a Binary files /dev/null and b/demos/plots/explicit/explicit Time_evolution_12 pi_100 steps_.png differ diff --git a/demos/plots/explicit/explicit Time_evolution_15 pi_100 steps_.png b/demos/plots/explicit/explicit Time_evolution_15 pi_100 steps_.png new file mode 100644 index 000000000..e50c89e6d Binary files /dev/null and b/demos/plots/explicit/explicit Time_evolution_15 pi_100 steps_.png differ diff --git a/demos/plots/explicit/explicit Time_evolution_4 pi_100 steps_.png b/demos/plots/explicit/explicit Time_evolution_4 pi_100 steps_.png new file mode 100644 index 000000000..76d8891f6 Binary files /dev/null and b/demos/plots/explicit/explicit Time_evolution_4 pi_100 steps_.png differ diff --git a/demos/plots/explicit/explicit Time_evolution_4 pi_1000 steps_.png b/demos/plots/explicit/explicit Time_evolution_4 pi_1000 steps_.png new file mode 100644 index 000000000..6ed0054bf Binary files /dev/null and b/demos/plots/explicit/explicit Time_evolution_4 pi_1000 steps_.png differ diff --git a/demos/plots/explicit/explicit Time_evolution_4 pi_10000 steps_.png b/demos/plots/explicit/explicit Time_evolution_4 pi_10000 steps_.png new file mode 100644 index 000000000..842077886 Binary files /dev/null and b/demos/plots/explicit/explicit Time_evolution_4 pi_10000 steps_.png differ diff --git a/demos/plots/explicit/explicit Time_evolution_4 pi_200 steps_.png b/demos/plots/explicit/explicit Time_evolution_4 pi_200 steps_.png new file mode 100644 index 000000000..c80bda326 Binary files /dev/null and b/demos/plots/explicit/explicit Time_evolution_4 pi_200 steps_.png differ diff --git a/demos/plots/explicit/explicit Time_evolution_60 pi_10000 steps_.png b/demos/plots/explicit/explicit Time_evolution_60 pi_10000 steps_.png new file mode 100644 index 000000000..bee511e64 Binary files /dev/null and b/demos/plots/explicit/explicit Time_evolution_60 pi_10000 steps_.png differ diff --git a/demos/plots/explicit/explicit Time_evolution_60 pi_1000000 steps_.png b/demos/plots/explicit/explicit Time_evolution_60 pi_1000000 steps_.png new file mode 100644 index 000000000..0d6559e0a Binary files /dev/null and b/demos/plots/explicit/explicit Time_evolution_60 pi_1000000 steps_.png differ diff --git a/demos/plots/explicit/explicit Time_evolution_60 pi_5000 steps_.png b/demos/plots/explicit/explicit Time_evolution_60 pi_5000 steps_.png new file mode 100644 index 000000000..10f3007d9 Binary files /dev/null and b/demos/plots/explicit/explicit Time_evolution_60 pi_5000 steps_.png differ diff --git a/demos/plots/implicit/.ipynb_checkpoints/implicit Phase_plot_4 pi_100 steps_-checkpoint.png b/demos/plots/implicit/.ipynb_checkpoints/implicit Phase_plot_4 pi_100 steps_-checkpoint.png new file mode 100644 index 000000000..f1ffe9433 Binary files /dev/null and b/demos/plots/implicit/.ipynb_checkpoints/implicit Phase_plot_4 pi_100 steps_-checkpoint.png differ diff --git a/demos/plots/implicit/implicit Phase_plot_4 pi_100 steps_.png b/demos/plots/implicit/implicit Phase_plot_4 pi_100 steps_.png new file mode 100644 index 000000000..f1ffe9433 Binary files /dev/null and b/demos/plots/implicit/implicit Phase_plot_4 pi_100 steps_.png differ diff --git a/demos/plots/implicit/implicit Time_evolution_4 pi_100 steps_.png b/demos/plots/implicit/implicit Time_evolution_4 pi_100 steps_.png new file mode 100644 index 000000000..797311564 Binary files /dev/null and b/demos/plots/implicit/implicit Time_evolution_4 pi_100 steps_.png differ diff --git a/demos/plots/improved/improved Phase_plot_10 pi_100 steps_.png b/demos/plots/improved/improved Phase_plot_10 pi_100 steps_.png new file mode 100644 index 000000000..46bd79994 Binary files /dev/null and b/demos/plots/improved/improved Phase_plot_10 pi_100 steps_.png differ diff --git a/demos/plots/improved/improved Phase_plot_4 pi_100 steps_.png b/demos/plots/improved/improved Phase_plot_4 pi_100 steps_.png new file mode 100644 index 000000000..967054e12 Binary files /dev/null and b/demos/plots/improved/improved Phase_plot_4 pi_100 steps_.png differ diff --git a/demos/plots/improved/improved Phase_plot_60 pi_1000 steps_.png b/demos/plots/improved/improved Phase_plot_60 pi_1000 steps_.png new file mode 100644 index 000000000..185a4d724 Binary files /dev/null and b/demos/plots/improved/improved Phase_plot_60 pi_1000 steps_.png differ diff --git a/demos/plots/improved/improved Phase_plot_60 pi_5000 steps_.png b/demos/plots/improved/improved Phase_plot_60 pi_5000 steps_.png new file mode 100644 index 000000000..02897931c Binary files /dev/null and b/demos/plots/improved/improved Phase_plot_60 pi_5000 steps_.png differ diff --git a/demos/plots/improved/improved Time_evolution_10 pi_100 steps_.png b/demos/plots/improved/improved Time_evolution_10 pi_100 steps_.png new file mode 100644 index 000000000..18693fa62 Binary files /dev/null and b/demos/plots/improved/improved Time_evolution_10 pi_100 steps_.png differ diff --git a/demos/plots/improved/improved Time_evolution_4 pi_100 steps_.png b/demos/plots/improved/improved Time_evolution_4 pi_100 steps_.png new file mode 100644 index 000000000..596a187ca Binary files /dev/null and b/demos/plots/improved/improved Time_evolution_4 pi_100 steps_.png differ diff --git a/demos/plots/improved/improved Time_evolution_60 pi_1000 steps_.png b/demos/plots/improved/improved Time_evolution_60 pi_1000 steps_.png new file mode 100644 index 000000000..db9906e04 Binary files /dev/null and b/demos/plots/improved/improved Time_evolution_60 pi_1000 steps_.png differ diff --git a/demos/plots/improved/improved Time_evolution_60 pi_5000 steps_.png b/demos/plots/improved/improved Time_evolution_60 pi_5000 steps_.png new file mode 100644 index 000000000..7cd78e73a Binary files /dev/null and b/demos/plots/improved/improved Time_evolution_60 pi_5000 steps_.png differ diff --git a/demos/test_autodiff.cpp b/demos/test_autodiff.cpp new file mode 100644 index 000000000..84f9a50c6 --- /dev/null +++ b/demos/test_autodiff.cpp @@ -0,0 +1,31 @@ +#include +#include "autodiff.hpp" + +using namespace ASC_ode; + +// f(x) = x^2 + 3 sin(x) +template +AutoDiff f(const AutoDiff& x) +{ + return pow(x, 2) + 3.0 * sin(x); +} + +int main() +{ + AutoDiff<1,double> x(1.0); + x.deriv()[0] = 1.0; // seed derivative + + auto y = f(x); + + std::cout << "f(1) = " << y.value() << "\n"; + std::cout << "f'(1) (AutoDiff) = " << y.deriv()[0] << "\n"; + + auto fnum = [](double z){ return z*z + 3.0*std::sin(z); }; + + double h = 1e-6; + double num = (fnum(1+h) - fnum(1-h)) / (2*h); + + std::cout << "f'(1) (numeric) = " << num << "\n"; + + return 0; +} diff --git a/demos/test_ode.cpp b/demos/test_ode.cpp index 97f635e40..cdf16fe56 100644 --- a/demos/test_ode.cpp +++ b/demos/test_ode.cpp @@ -1,60 +1,100 @@ #include -#include - -#include -#include +#include +#include +#include "nonlinfunc.hpp" +#include "timestepper.hpp" using namespace ASC_ode; - +// ======================================================== +// Mass–Spring System: y' = f(y) +// y = [x, v] +// x' = v +// v' = -(k/m) * x +// ======================================================== class MassSpring : public NonlinearFunction { private: - double mass; - double stiffness; + double mass; + double stiffness; public: - MassSpring(double m, double k) : mass(m), stiffness(k) {} - - size_t dimX() const override { return 2; } - size_t dimF() const override { return 2; } - - void evaluate (VectorView x, VectorView f) const override - { - f(0) = x(1); - f(1) = -stiffness/mass*x(0); - } - - void evaluateDeriv (VectorView x, MatrixView df) const override - { - df = 0.0; - df(0,1) = 1; - df(1,0) = -stiffness/mass; - } + MassSpring(double m, double k) : mass(m), stiffness(k) {} + + size_t dimX() const override { return 2; } + size_t dimF() const override { return 2; } + + void evaluate(VectorView x, VectorView f) const override + { + f(0) = x(1); + f(1) = -stiffness / mass * x(0); + } + + void evaluateDeriv(VectorView x, MatrixView df) const override + { + df = 0.0; + df(0, 1) = 1.0; + df(1, 0) = -stiffness / mass; + } }; -int main() +int main(int argc, char* argv[]) { - double tend = 4*M_PI; - int steps = 100; - double tau = tend/steps; - - Vector<> y = { 1, 0 }; // initializer list - auto rhs = std::make_shared(1.0, 1.0); - - ExplicitEuler stepper(rhs); - // ImplicitEuler stepper(rhs); - - std::ofstream outfile ("output_test_ode.txt"); - std::cout << 0.0 << " " << y(0) << " " << y(1) << std::endl; - outfile << 0.0 << " " << y(0) << " " << y(1) << std::endl; - - for (int i = 0; i < steps; i++) - { - stepper.DoStep(tau, y); - - std::cout << (i+1) * tau << " " << y(0) << " " << y(1) << std::endl; - outfile << (i+1) * tau << " " << y(0) << " " << y(1) << std::endl; - } + if (argc < 4) + { + std::cout << "Usage: ./test_ode T_relative steps method\n"; + std::cout << "Example: ./test_ode 4 100 RK2\n"; + return 1; + } + + double tend_relative = atof(argv[1]); + int steps = atoi(argv[2]); + std::string algorithm = argv[3]; + + double tend = tend_relative * M_PI; + double tau = tend / steps; + + Vector<> y = {1.0, 0.0}; + auto rhs = std::make_shared(1.0, 1.0); + + // ======================================================== + // Choose timestepper + // ======================================================== + std::unique_ptr stepper; + + if (algorithm == "explicit") + stepper = std::make_unique(rhs); + else if (algorithm == "improved") + stepper = std::make_unique(rhs); + else if (algorithm == "implicit") + stepper = std::make_unique(rhs); + else if (algorithm == "CN") + stepper = std::make_unique(rhs); + else if (algorithm == "RK2") + stepper = std::make_unique(rhs); + else + { + std::cout << "Choose method: explicit / improved / implicit / CN / RK2\n"; + return 1; + } + + // ======================================================== + // Write output file + // ======================================================== + std::ofstream outfile("output_test_ode.txt"); + + std::cout << 0.0 << " " << y(0) << " " << y(1) << "\n"; + outfile << 0.0 << " " << y(0) << " " << y(1) << "\n"; + + for (int i = 0; i < steps; i++) + { + stepper->DoStep(tau, y); + + double t = (i + 1) * tau; + std::cout << t << " " << y(0) << " " << y(1) << "\n"; + outfile << t << " " << y(0) << " " << y(1) << "\n"; + } + + return 0; } diff --git a/demos/test_pendulum_autodiff.cpp b/demos/test_pendulum_autodiff.cpp new file mode 100644 index 000000000..c3911038e --- /dev/null +++ b/demos/test_pendulum_autodiff.cpp @@ -0,0 +1,34 @@ +#include +#include +#include "../src/autodiff.hpp" + +using namespace ASC_ode; + +// Pendulum equation: f(theta) = - (g/L) * sin(theta) +constexpr double g = 9.81; +constexpr double L = 1.0; + +int main() +{ + using AD = AutoDiff<1,double>; + + // AutoDiff variable: theta = 0.5 rad + AD theta(0.5); + theta.deriv()[0] = 1.0; // d(theta)/d(theta) = 1 + + // Pendulum function + AD f = -(g/L) * sin(theta); + + // Numerical derivative for comparison + double h = 1e-6; + double f_num = -(g/L) * sin(0.5); + double f_plus = -(g/L) * sin(0.5 + h); + double df_numeric = (f_plus - f_num) / h; + + std::cout << "Pendulum AutoDiff test\n"; + std::cout << "Value f(theta) = " << f.value() << "\n"; + std::cout << "AutoDiff derivative df/dtheta = " << f.deriv()[0] << "\n"; + std::cout << "Numeric derivative df/dtheta = " << df_numeric << "\n"; + + return 0; +} diff --git a/demos/test_rc.cpp b/demos/test_rc.cpp new file mode 100644 index 000000000..70e290dea --- /dev/null +++ b/demos/test_rc.cpp @@ -0,0 +1,63 @@ +#include +#include +#include +#include +#include "timestepper.hpp" +#include "rccircuit.hpp" + +using namespace ASC_ode; + +int main(int argc, char* argv[]) +{ + if(argc < 4){ + std::cout << "Usage: ./test_rc tend steps method\n"; + return 1; + } + + double tend = atof(argv[1]); + int steps = atoi(argv[2]); + std::string method = argv[3]; + + double tau = tend / steps; + + // initial state: UC(0) = 0, t(0) = 0 + Vector<> y = {0.0, 0.0}; + + // RC circuit: R = 100 Ohm, C = 1e-6 F + auto rhs = std::make_shared(100.0, 1e-6); + + std::unique_ptr stepper; + + if(method == "explicit") + stepper = std::make_unique(rhs); + else if(method == "improved") + stepper = std::make_unique(rhs); + else if(method == "implicit") + stepper = std::make_unique(rhs); + else if(method == "CN") + stepper = std::make_unique(rhs); + else { + std::cout << "Choose: explicit / improved / implicit / CN\n"; + return 1; + } + + // ----------------------------- + // Output file per method: + // rc_explicit.txt + // rc_improved.txt + // rc_implicit.txt + // rc_CN.txt + // ----------------------------- + std::string fname = "rc_" + method + ".txt"; + std::ofstream out(fname); + + out << 0 << " " << y(0) << " " << y(1) << "\n"; + + for(int i = 0; i < steps; i++){ + stepper->DoStep(tau, y); + out << (i+1) * tau << " " << y(0) << " " << y(1) << "\n"; + } + + std::cout << "Saved: " << fname << "\n"; + return 0; +} diff --git a/nanoblas b/nanoblas index 8204aef71..f935be2ef 160000 --- a/nanoblas +++ b/nanoblas @@ -1 +1 @@ -Subproject commit 8204aef7160cee3e9c18265f96b7211e7c479ad2 +Subproject commit f935be2ef996f8693c49e1203014355615d762fb diff --git a/report_images/Phase_plot_standard.png b/report_images/Phase_plot_standard.png new file mode 100644 index 000000000..3d1db1982 Binary files /dev/null and b/report_images/Phase_plot_standard.png differ diff --git a/report_images/RC_Time_evolution_tend_0.01_steps_2000.png b/report_images/RC_Time_evolution_tend_0.01_steps_2000.png new file mode 100644 index 000000000..49b3de8f5 Binary files /dev/null and b/report_images/RC_Time_evolution_tend_0.01_steps_2000.png differ diff --git a/report_images/RC_Time_evolution_zoom_tend_0.01_steps_2000.png b/report_images/RC_Time_evolution_zoom_tend_0.01_steps_2000.png new file mode 100644 index 000000000..af9e12917 Binary files /dev/null and b/report_images/RC_Time_evolution_zoom_tend_0.01_steps_2000.png differ diff --git a/report_images/Time_evolution_standard.png b/report_images/Time_evolution_standard.png new file mode 100644 index 000000000..f9bcc6824 Binary files /dev/null and b/report_images/Time_evolution_standard.png differ diff --git a/runmassspring.sh b/runmassspring.sh new file mode 100755 index 000000000..a42139b39 --- /dev/null +++ b/runmassspring.sh @@ -0,0 +1,13 @@ +#!/bin/bash + +tend_relative=$1 #standard: 4 +steps=$2 #standard: 100 +algorithm=$3 + +cd build +./test_ode > ../demos/output_test_ode.txt $tend_relative $steps $algorithm +cd ../demos +python3 plotmassspring.py --tend_relative $1 --steps $2 --algorithm $3 +mv "Phase_plot.png" "plots/$3/$3 Phase_plot_$1 pi_$2 steps_.png" +mv "Time_evolution.png" "plots/$3/$3 Time_evolution_$1 pi_$2 steps_.png" +cd .. \ No newline at end of file diff --git a/src/autodiff.hpp b/src/autodiff.hpp index 40a46e7e5..3ca6a097f 100644 --- a/src/autodiff.hpp +++ b/src/autodiff.hpp @@ -1,110 +1,251 @@ #ifndef AUTODIFF_HPP #define AUTODIFF_HPP -#include -#include -#include -#include - +#include +#include +#include namespace ASC_ode { - template - class Variable - { - private: - T m_val; - public: - Variable (T v) : m_val(v) {} - T value() const { return m_val; } - }; +// ===================================================== +// Variable class (template for independent variable xi) +// ===================================================== +template +class Variable +{ +private: + T m_val; + +public: + Variable(T v) : m_val(v) {} + T value() const { return m_val; } +}; - template - auto derivative (T v, size_t /*index*/) { return T(0); } +// fallback for normal double +template +T derivative(T, size_t) { return T(0); } - template - class AutoDiff - { - private: +// ===================================================== +// AutoDiff class +// ===================================================== +template +class AutoDiff +{ +private: T m_val; std::array m_deriv; - public: - AutoDiff () : m_val(0), m_deriv{} {} - AutoDiff (T v) : m_val(v), m_deriv{} + +public: + AutoDiff() : m_val(0.0), m_deriv{} {} + + AutoDiff(T v) : m_val(v), m_deriv{} { - for (size_t i = 0; i < N; i++) - m_deriv[i] = derivative(v, i); + for (size_t i = 0; i < N; i++) + m_deriv[i] = derivative(v, i); } - + template - AutoDiff (Variable var) : m_val(var.value()), m_deriv{} + AutoDiff(Variable var) : m_val(var.value()), m_deriv{} { - m_deriv[I] = 1.0; + m_deriv[I] = 1.0; } T value() const { return m_val; } - std::array& deriv() { return m_deriv; } - const std::array& deriv() const { return m_deriv; } - }; + std::array &deriv() { return m_deriv; } + const std::array &deriv() const { return m_deriv; } +}; - template - auto derivative (AutoDiff v, size_t index) - { +// ===================================================== +// derivative() for AutoDiff type +// ===================================================== +template +auto derivative(const AutoDiff &v, size_t index) +{ return v.deriv()[index]; - } - +} - template - std::ostream & operator<< (std::ostream& os, const AutoDiff& ad) - { +// ===================================================== +// Printing operator +// ===================================================== +template +std::ostream &operator<<(std::ostream &os, const AutoDiff &ad) +{ os << "Value: " << ad.value() << ", Deriv: ["; for (size_t i = 0; i < N; i++) { - os << ad.deriv()[i]; - if (i < N - 1) os << ", "; + os << ad.deriv()[i]; + if (i < N - 1) os << ", "; } os << "]"; return os; - } - - template - AutoDiff operator+ (const AutoDiff& a, const AutoDiff& b) - { - AutoDiff result(a.value() + b.value()); - for (size_t i = 0; i < N; i++) - result.deriv()[i] = a.deriv()[i] + b.deriv()[i]; - return result; - } - - template - auto operator+ (T a, const AutoDiff& b) { return AutoDiff(a) + b; } - - - template - AutoDiff operator* (const AutoDiff& a, const AutoDiff& b) - { - AutoDiff result(a.value() * b.value()); - for (size_t i = 0; i < N; i++) - result.deriv()[i] = a.deriv()[i] * b.value() + a.value() * b.deriv()[i]; - return result; - } - - using std::sin; - using std::cos; - - template - AutoDiff sin(const AutoDiff &a) - { - AutoDiff result(sin(a.value())); - for (size_t i = 0; i < N; i++) - result.deriv()[i] = cos(a.value()) * a.deriv()[i]; - return result; - } +} + + +// ===================================================== +// Addition +// ===================================================== +template +AutoDiff operator+(const AutoDiff &a, const AutoDiff &b) +{ + AutoDiff r(a.value() + b.value()); + for (size_t i = 0; i < N; i++) + r.deriv()[i] = a.deriv()[i] + b.deriv()[i]; + return r; +} + +template +AutoDiff operator+(T a, const AutoDiff &b) +{ + return AutoDiff(a) + b; +} + +template +AutoDiff operator+(const AutoDiff &a, T b) +{ + return a + AutoDiff(b); +} + + +// ===================================================== +// Subtraction +// ===================================================== +template +AutoDiff operator-(const AutoDiff &a) +{ + AutoDiff r(-a.value()); + for (size_t i = 0; i < N; i++) + r.deriv()[i] = -a.deriv()[i]; + return r; +} + +template +AutoDiff operator-(const AutoDiff &a, const AutoDiff &b) +{ + AutoDiff r(a.value() - b.value()); + for (size_t i = 0; i < N; i++) + r.deriv()[i] = a.deriv()[i] - b.deriv()[i]; + return r; +} + + +// ===================================================== +// Multiplication (3 overloads) +// ===================================================== +template +AutoDiff operator*(const AutoDiff &a, const AutoDiff &b) +{ + AutoDiff r(a.value() * b.value()); + for (size_t i = 0; i < N; i++) + r.deriv()[i] = a.deriv()[i] * b.value() + a.value() * b.deriv()[i]; + return r; +} + +template +AutoDiff operator*(T a, const AutoDiff &b) +{ + return AutoDiff(a) * b; +} +template +AutoDiff operator*(const AutoDiff &a, T b) +{ + return a * AutoDiff(b); +} + + +// ===================================================== +// Division (3 full overloads) +// ===================================================== +template +AutoDiff operator/(const AutoDiff &a, const AutoDiff &b) +{ + AutoDiff r(a.value() / b.value()); + T b2 = b.value() * b.value(); + + for (size_t i = 0; i < N; i++) + r.deriv()[i] = (a.deriv()[i] * b.value() - a.value() * b.deriv()[i]) / b2; + + return r; +} + +template +AutoDiff operator/(const AutoDiff &a, T c) +{ + AutoDiff r(a.value() / c); + for (size_t i = 0; i < N; i++) + r.deriv()[i] = a.deriv()[i] / c; + return r; +} + +template +AutoDiff operator/(T c, const AutoDiff &a) +{ + AutoDiff r(c / a.value()); + T a2 = a.value() * a.value(); + + for (size_t i = 0; i < N; i++) + r.deriv()[i] = -c * a.deriv()[i] / a2; + + return r; +} + + +// ===================================================== +// Elementary functions +// ===================================================== +using std::cos; +using std::exp; +using std::log; +using std::sin; + +template +AutoDiff sin(const AutoDiff &a) +{ + AutoDiff r(sin(a.value())); + for (size_t i = 0; i < N; i++) + r.deriv()[i] = cos(a.value()) * a.deriv()[i]; + return r; +} + +template +AutoDiff cos(const AutoDiff &a) +{ + AutoDiff r(cos(a.value())); + for (size_t i = 0; i < N; i++) + r.deriv()[i] = -sin(a.value()) * a.deriv()[i]; + return r; +} + +template +AutoDiff exp(const AutoDiff &a) +{ + AutoDiff r(exp(a.value())); + for (size_t i = 0; i < N; i++) + r.deriv()[i] = exp(a.value()) * a.deriv()[i]; + return r; +} + +template +AutoDiff log(const AutoDiff &a) +{ + AutoDiff r(log(a.value())); + for (size_t i = 0; i < N; i++) + r.deriv()[i] = a.deriv()[i] / a.value(); + return r; +} + +template +AutoDiff pow(const AutoDiff &a, double p) +{ + AutoDiff r(std::pow(a.value(), p)); + for (size_t i = 0; i < N; i++) + r.deriv()[i] = p * std::pow(a.value(), p - 1) * a.deriv()[i]; + return r; +} } // namespace ASC_ode diff --git a/src/implicitRK.hpp b/src/implicitRK.hpp new file mode 100644 index 000000000..e7434b334 --- /dev/null +++ b/src/implicitRK.hpp @@ -0,0 +1,231 @@ +#ifndef IMPLICITRK_HPP +#define IMPLICITRK_HPP + +#include +#include +#include + +namespace ASC_ode { + using namespace nanoblas; + + + + class ImplicitRungeKutta : public TimeStepper + { + Matrix<> m_a; + Vector<> m_b, m_c; + std::shared_ptr m_equ; + std::shared_ptr m_tau; + std::shared_ptr m_yold; + int m_stages; + int m_n; + Vector<> m_k, m_y; + public: + ImplicitRungeKutta(std::shared_ptr rhs, + const Matrix<> &a, const Vector<> &b, const Vector<> &c) + : TimeStepper(rhs), m_a(a), m_b(b), m_c(c), + m_tau(std::make_shared(0.0)), + m_stages(c.size()), m_n(rhs->dimX()), m_k(m_stages*m_n), m_y(m_stages*m_n) + { + auto multiple_rhs = make_shared(rhs, m_stages); + m_yold = std::make_shared(m_stages*m_n); + auto knew = std::make_shared(m_stages*m_n); + m_equ = knew - Compose(multiple_rhs, m_yold+m_tau*std::make_shared(a, m_n)); + } + + void DoStep(double tau, VectorView y) override + { + for (int j = 0; j < m_stages; j++) + m_y.range(j*m_n, (j+1)*m_n) = y; + m_yold->set(m_y); + + m_tau->set(tau); + m_k = 0.0; + NewtonSolver(m_equ, m_k); + + for (int j = 0; j < m_stages; j++) + y += tau * m_b(j) * m_k.range(j*m_n, (j+1)*m_n); + } + }; + + + + + + + +Matrix Gauss2a { { 0.25, 0.25 - sqrt(3)/6 }, { 0.25 + sqrt(3)/6, 0.25 } }; +Vector<> Gauss2b { 0.5, 0.5 }; +Vector<> Gauss2c { 0.5 - sqrt(3)/6, 0.5 + sqrt(3)/6 }; + + +Vector<> Gauss3c { 0.5 - sqrt(15)/10, 0.5, 0.5+sqrt(15)/10 }; + + +// codes from Numerical Recipes, https://numerical.recipes/book.html + +// Gauss integration on [0,1] +void GaussLegendre(VectorView<> x, VectorView<> w) +// Given the lower and upper limits of integration x1 and x2, this routine returns arrays x[0..n-1] +// and w[0..n-1] of length n, containing the abscissas and weights of the Gauss-Legendre n-point +// quadrature formula. + { + double x1 = 0; + double x2 = 1; + const double EPS=1.0e-14; // EPS is the relative precision. + double z1,z,xm,xl,pp,p3,p2,p1; + int n=x.size(); + int m=(n+1)/2; // The roots are symmetric in the interval, so + xm=0.5*(x2+x1); // we only have to find half of them. + xl=0.5*(x2-x1); + for (int i=0;i EPS); + x[i]=xm-xl*z; // Scale the root to the desired interval, + x[n-1-i]=xm+xl*z; // and put in its symmetric counterpart. + w[i]=2.0*xl/((1.0-z*z)*pp*pp); // Compute the weight + w[n-1-i]=w[i]; // and its symmetric counterpart. + } + } + + +void GaussJacobi (VectorView<> x, VectorView<> w, const double alf, const double bet) +// Given alf and bet, the parameters ˛ and ˇ of the Jacobi polynomials, this routine returns +// arrays x[0..n-1] and w[0..n-1] containing the abscissas and weights of the n-point GaussJacobi quadrature formula. The largest abscissa is returned in x[0], the smallest in x[n-1]. +{ + const int MAXIT=10; + const double EPS=1.0e-14; // EPS is the relative precision. + int i,its,j; + double alfbet,an,bn,r1,r2,r3; + double a,b,c,p1,p2,p3,pp,temp,z,z1; + int n=x.size(); + for (i=0;i MAXIT) throw("too many iterations in gaujac"); + x[i]=z; // Store the root and the weight. + w[i]=exp(std::lgamma(alf+n)+std::lgamma(bet+n)-std::lgamma(n+1.0)- + std::lgamma(n+alfbet+1.0))*temp*pow(2.0,alfbet)/(pp*p2); + } +} + + + +/* + given Runge-Kutta nodes c, compute the coefficients a and b +*/ +auto ComputeABfromC (const Vector<> & c) +{ + int s = c.size(); + Matrix<> M(s, s); + Vector<> tmp(s); + + for (int i = 0; i < s; i++) + for (int j = 0; j < s; j++) + M(i,j) = std::pow(c(j), i); + + calcInverse(M); + // M = LapackLU(M).inverse(); + + for (int i = 0; i < s; i++) + tmp(i) = 1.0 / (i+1); + + Vector<> b = M * tmp; + Matrix a(s,s); + + for (int j = 0; j < s; j++) + { + for (int i = 0; i < s; i++) + tmp(i) = std::pow(c(j),i+1) / (i+1); + a.row(j) = M * tmp; + } + /* + std::cout << "b = " << b << std::endl; + std::cout << "a = " << a << std::endl; + */ + return std::tuple { a, b }; +} + + +void GaussRadau (VectorView<> x, VectorView<> w) +{ + GaussJacobi (x.range(0, x.size()-1), + w.range(0, w.size()-1), 1, 0); + for (int i = 0; i < x.size()-1; i++) + { + x(i) = 0.5*(x(i)+1); + w(i) *= 0.5; + } + x(x.size()-1) = 1.0; + double sum = 0; + for (int i = 0; i < x.size()-1; i++) + sum += w(i); + w(x.size()-1) = 1-sum; +} +} + +#endif // IMPLICITRK_HPP \ No newline at end of file diff --git a/src/massspring.cpp b/src/massspring.cpp new file mode 100644 index 000000000..1d8314891 --- /dev/null +++ b/src/massspring.cpp @@ -0,0 +1,5 @@ +#include + +int main(){ + +} diff --git a/src/nonlinfunc.hpp b/src/nonlinfunc.hpp index 76a19957c..3f69000b0 100644 --- a/src/nonlinfunc.hpp +++ b/src/nonlinfunc.hpp @@ -244,6 +244,61 @@ namespace ASC_ode }; + class MultipleFunc : public NonlinearFunction + { + std::shared_ptr func; + size_t num, fdimx, fdimf; + public: + MultipleFunc (std::shared_ptr _func, int _num) + : func(_func), num(_num) + { + fdimx = func->dimX(); + fdimf = func->dimF(); + } + + virtual size_t dimX() const override { return num * fdimx; } + virtual size_t dimF() const override{ return num * fdimf; } + virtual void evaluate (VectorView x, VectorView f) const override + { + for (size_t i = 0; i < num; i++) + func->evaluate(x.range(i*fdimx, (i+1)*fdimx), + f.range(i*fdimf, (i+1)*fdimf)); + } + virtual void evaluateDeriv (VectorView x, MatrixView df) const override + { + df = 0.0; + for (size_t i = 0; i < num; i++) + func->evaluateDeriv(x.range(i*fdimx, (i+1)*fdimx), + df.rows(i*fdimf, (i+1)*fdimf).cols(i*fdimx, (i+1)*fdimx)); + } + }; + + + class MatVecFunc : public NonlinearFunction + { + Matrix<> m_a; + size_t m_n; + public: + MatVecFunc (Matrix<> a, size_t n) + : m_a(a), m_n(n) { } + + virtual size_t dimX() const override { return m_n*m_a.rows(); } + virtual size_t dimF() const override { return m_n*m_a.cols(); } + virtual void evaluate (VectorView x, VectorView f) const override + { + MatrixView mx(m_a.cols(), m_n, m_n, x.data()); + MatrixView mf(m_a.rows(), m_n, m_n, f.data()); + mf = m_a * mx; + } + virtual void evaluateDeriv (VectorView x, MatrixView df) const override + { + df = 0.0; + for (size_t i = 0; i < m_a.rows(); i++) + for (size_t j = 0; j < m_a.cols(); j++) + df.rows(i*m_n, (i+1)*m_n).cols(j*m_n, (j+1)*m_n).diag() = m_a(i,j); + } + }; + } #endif diff --git a/src/rccircuit.hpp b/src/rccircuit.hpp new file mode 100644 index 000000000..f1a6d7e9c --- /dev/null +++ b/src/rccircuit.hpp @@ -0,0 +1,35 @@ +#pragma once +#include "nonlinfunc.hpp" +#include + +namespace ASC_ode { + +class RCCircuit : public NonlinearFunction +{ + double R, C; + +public: + RCCircuit(double R_, double C_) : R(R_), C(C_) {} + + size_t dimX() const override { return 2; } + size_t dimF() const override { return 2; } + + // x(0) = U_C (capacitor voltage) + // x(1) = t (time variable) + void evaluate(VectorView x, VectorView f) const override { + double UC = x(0); + double t = x(1); + double U0 = std::cos(100.0 * M_PI * t); + + f(0) = (U0 - UC) / (R*C); // dUC/dt + f(1) = 1.0; // dt/dt = 1 + } + + void evaluateDeriv(VectorView x, MatrixView df) const override { + df = 0.0; + df(0,0) = -1.0/(R*C); + df(0,1) = (-100.0*M_PI * std::sin(100.0 * M_PI * x(1))) / (R*C); + } +}; + +} diff --git a/src/timestepper.hpp b/src/timestepper.hpp index c7802243b..bc97804e4 100644 --- a/src/timestepper.hpp +++ b/src/timestepper.hpp @@ -33,6 +33,21 @@ namespace ASC_ode } }; + class ImprovedEuler : public TimeStepper + { + Vector<> m_vecf; + public: + ImprovedEuler(std::shared_ptr rhs) + : TimeStepper(rhs), m_vecf(rhs->dimF()) {} + void DoStep(double tau, VectorView y) override + { + Vector y_tilde(y.size()); + y_tilde = y + 0.5 *tau * m_vecf; + this->m_rhs->evaluate(y_tilde, m_vecf); + y += tau * m_vecf; + } + }; + class ImplicitEuler : public TimeStepper { std::shared_ptr m_equ; @@ -54,11 +69,104 @@ namespace ASC_ode NewtonSolver(m_equ, y); } }; + class RungeKutta2 : public TimeStepper + { + public: + using TimeStepper::TimeStepper; + + void DoStep(double tau, VectorView y) override + { + const size_t n = m_rhs->dimX(); + + Vector k1(n); + Vector k2(n); + Vector ytmp(n); + + // k1 = f(y_n) + m_rhs->evaluate(y, k1); + + // ytmp = y_n + (tau/2) * k1 + for (size_t i = 0; i < n; ++i) + ytmp(i) = y(i) + 0.5 * tau * k1(i); + + // k2 = f(y_n + tau/2 * k1) + m_rhs->evaluate(ytmp, k2); + + // y_{n+1} = y_n + tau * k2 + for (size_t i = 0; i < n; ++i) + y(i) += tau * k2(i); + } + }; - + class RungeKutta4 : public TimeStepper + { + public: + using TimeStepper::TimeStepper; + + void DoStep(double tau, VectorView y) override + { + const size_t n = m_rhs->dimX(); + Vector k1(n); + Vector k2(n); + Vector k3(n); + Vector k4(n); + Vector ytmp(n); + + // k1 = f(y_n) + m_rhs->evaluate(y, k1); + + // k2 = f(y_n + tau/2 * k1) + for (size_t i = 0; i < n; ++i) + ytmp(i) = y(i) + 0.5 * tau * k1(i); + m_rhs->evaluate(ytmp, k2); + + // k3 = f(y_n + tau/2 * k2) + for (size_t i = 0; i < n; ++i) + ytmp(i) = y(i) + 0.5 * tau * k2(i); + m_rhs->evaluate(ytmp, k3); + + // k4 = f(y_n + tau * k3) + for (size_t i = 0; i < n; ++i) + ytmp(i) = y(i) + tau * k3(i); + m_rhs->evaluate(ytmp, k4); + + // y_{n+1} = y_n + tau/6 * (k1 + 2k2 + 2k3 + k4) + for (size_t i = 0; i < n; ++i) + { + y(i) += (tau / 6.0) * + (k1(i) + 2.0 * k2(i) + 2.0 * k3(i) + k4(i)); + } + } + }; + + class CrankNicolson : public TimeStepper + { + std::shared_ptr m_equ; + std::shared_ptr m_tau; + std::shared_ptr m_yold; + public: + CrankNicolson(std::shared_ptr rhs) + : TimeStepper(rhs), m_tau(std::make_shared(0.0)) + { + m_yold = std::make_shared(rhs->dimX()); + auto ynew = std::make_shared(rhs->dimX()); + auto f_old = std::make_shared(rhs, m_yold); + + m_equ = ynew - m_yold - m_tau * (f_old + m_rhs); + } + + void DoStep(double tau, VectorView y) override + { + m_yold->set(y); + m_tau->set(0.5 * tau); + NewtonSolver(m_equ, y); + } + + }; } + #endif