Add HW17.4, HW18, HW19 implementations + full report

.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 ..

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)

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

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.
+<h4 style="text-align:center;">Time Evolution (Mass–Spring)</h4>
+<p align="center">
+  <img src="report_images/Time_evolution_standard.png" width="450">
+</p>
+
+<h4 style="text-align:center;">Phase Plot (Mass–Spring)</h4>
+<p align="center">
+  <img src="report_images/Phase_plot_standard.png" width="450">
+</p>
+
+#### **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**
+
+<h4 style="text-align:center;">RC Circuit – Full Time Evolution</h4>
+<p align="center">
+  <img src="report_images/RC_Time_evolution_tend_0.01_steps_2000.png" width="500">
+</p>
+
+
+<h4 style="text-align:center;">RC Circuit – Zoomed View (Stability)</h4>
+<p align="center">
+  <img src="report_images/RC_Time_evolution_zoom_tend_0.01_steps_2000.png" width="500">
+</p>
+
+#### **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<N,T>`**, 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<I>`
+- 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  
+
+
+
+---

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:
+
+<div style="display: inline-block; text-align: center; margin-right: 10px;">
+  <img src="demos/plots/explicit/explicit Time_evolution_4 pi_100 steps_.png" alt="Time_evolution" width="300">
+</div>
+<div style="display: inline-block; text-align: center;">
+  <img src="demos/plots/explicit/explicit Phase_plot_4 pi_100 steps_.png" alt="Phase_plot" width="300">
+</div>
+
+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:
+
+<div style="display: inline-block; text-align: center; margin-right: 10px;">
+  <img src="demos/plots/explicit/explicit Time_evolution_4 pi_10000 steps_.png" alt="Time_evolution" width="300">
+</div>
+<div style="display: inline-block; text-align: center;">
+  <img src="demos/plots/explicit/explicit Phase_plot_4 pi_10000 steps_.png" alt="Phase_plot" width="300">
+</div>
+
+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.
+
+<div style="display: inline-block; text-align: center; margin-right: 10px;">
+  <img src="demos/plots/explicit/explicit Time_evolution_60 pi_10000 steps_.png" alt="Time_evolution" width="300">
+</div>
+<div style="display: inline-block; text-align: center;">
+  <img src="demos/plots/explicit/explicit Phase_plot_60 pi_10000 steps_.png" alt="Phase_plot" width="300">
+</div>
+
+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}}
+$$
+
+<div style="display: inline-block; text-align: center; margin-right: 10px;">
+  <img src="demos/plots/explicit/explicit Time_evolution_60 pi_1000000 steps_.png" alt="Time_evolution" width="300">
+</div>
+<div style="display: inline-block; text-align: center;">
+  <img src="demos/plots/explicit/explicit Phase_plot_60 pi_1000000 steps_.png" alt="Phase_plot" width="300">
+</div>
+
+
+### 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.
+
+
+
+<div style="display: inline-block; text-align: center; margin-right: 10px;">
+  <img src="demos/plots/explicit/explicit Phase_plot_10 pi_100 steps_.png" alt="Time_evolution" width="300">
+</div>
+<div style="display: inline-block; text-align: center;">
+  <img src="demos/plots/improved/improved Phase_plot_10 pi_100 steps_.png" alt="Phase_plot" width="300">  
+</div>
+
+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.
+
+<div style="display: inline-block; text-align: center; margin-right: 10px;">
+  <img src="demos/plots/explicit/explicit Time_evolution_60 pi_5000 steps_.png" alt="Time_evolution" width="300">
+</div>
+<div style="display: inline-block; text-align: center;">
+  <img src="demos/plots/improved/improved Time_evolution_60 pi_5000 steps_.png" alt="Phase_plot" width="300">
+</div>
+
+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
+$$
+
+
+<div style="display: inline-block; text-align: center; margin-right: 10px;">
+  <img src="demos/plots/explicit/explicit Phase_plot_60 pi_5000 steps_.png" alt="Time_evolution" width="300">
+</div>
+<div style="display: inline-block; text-align: center;">
+  <img src="demos/plots/improved/improved Phase_plot_60 pi_5000 steps_.png" alt="Phase_plot" width="300">
+</div>