Skip to content

Commit d202089

Browse files
tlunettlunet
committed
TL: almost finished chap2
1 parent 6088a31 commit d202089

33 files changed

+737
-207
lines changed

‎.gitignore

+2-1
Original file line numberDiff line numberDiff line change
@@ -2,4 +2,5 @@
22
*.png
33
*.pyc
44
*.m~
5-
.spyproject
5+
.spyproject
6+
*.eps

‎README.md

+4-4
Original file line numberDiff line numberDiff line change
@@ -28,11 +28,11 @@ The folder organization follows the chapter structure of the book, as follow :
2828
- [1.4 : Historical Overview](./chap1/README.md#section-14--historical-overview)
2929
- [1.5 : Problems](./chap1/sec1.5/README.md)
3030

31-
## Chapter 2 : Multiple Shooting Type Methods
31+
## [Chapter 2 : Multiple Shooting Type Methods](./chap2/README.md)
3232

33-
- [2.1 : Idea of Nievergelt in 1964](./chap2/sec2.1/README.md)
34-
- [2.2 : Multiple Shooting Methods in Time](./chap2/sec2.2/README.md)
35-
- [2.3 : The Parareal Algorithm](./chap2/sec2.3/README.md)
33+
- [2.1 : Idea of Nievergelt in 1964](./chap2/README.md#section-21--idea-of-nievergelt-in-1964)
34+
- [2.2 : Multiple Shooting Methods in Time](./chap2/README.md#section-22--multiple-shooting-methods-in-time)
35+
- [2.3 : The Parareal Algorithm](./chap2/README.md#section-23--the-parareal-algorithm)
3636
- [2.4 : Problems](./chap2/sec2.4/README.md)
3737

3838
## Chapter 3 : Waveform Relaxation and Domain Decomposition

‎chap1/README.md

+1-1
Original file line numberDiff line numberDiff line change
@@ -23,4 +23,4 @@
2323

2424
## [Section 1.5 : Problems](sec1.5/README.md)
2525

26-
[📖 Book Content](../README.md)
26+
[📖 Table of Content](../README.md)

‎chap1/sec1.5/README.md

+1-1
Original file line numberDiff line numberDiff line change
@@ -11,4 +11,4 @@
1111
- [Problem 1.9 : Implementation and testing of a wave equation solver](./prob1.9/README.md)
1212
- [Problem 1.10 : Analysis of the Newton's method](./prob1.10/README.md)
1313

14-
[📖 Book Content](../../README.md) | [📄 Chapter Content](../README.md)
14+
[📖 Table of Content](../../README.md) | [📄 Chapter Sections](../README.md)

‎chap2/README.md

+34
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
# Chapter 2 : Multiple Shooting Type Methods
2+
3+
## Section 2.1 : Idea of Nievergelt in 1964
4+
5+
| 📜 Description | 🧮 Matlab | 🐍 Python |
6+
| :--- | :--- | :--- |
7+
| implementation of the Nievergelt method | [Nievergelt.m](./sec2.1/Nievergelt.m) | [Nievergelt.py](./sec2.1/Nievergelt.py) |
8+
| generic implementation of Forward Euler | [ForwardEuler.m](./sec2.1/ForwardEuler.m), [UForwardEuler.m](./sec2.1/UForwardEuler.m), [exampleODELinear.m](sec2.1/exampleODELinear.m) | [ForwardEuler.py](./sec2.1/ForwardEuler.py), [exampleODELinear.py](sec2.1/exampleODELinear.py) |
9+
| Utility function to find the index of the closest interval | [FindClosestInterval.m](./sec2.1/FindClosestInterval.m) | _in `Nievergelt.py`_ |
10+
| Example code for Nievergelt method with the linear problem | [exampleNievergelt.m](./sec2.1/exampleNievergelt.m) | [exampleNievergelt.py](./sec2.1/exampleNievergelt.py) |
11+
| Script to generate the figure for the linear problem | [exampleNievergeltLinearFig.m](./sec2.1/exampleNievergeltLinearFig.m) | [exampleNievergeltLinearFig.py](./sec2.1/exampleNievergeltLinearFig.py) |
12+
| Script to generate the figure for the non-linear prolem | [exampleNievergeltNonLinearFig.m](./sec2.1/exampleNievergeltNonLinearFig.m) | [exampleNievergeltNonLinearFig.py](./sec2.1/exampleNievergeltNonLinearFig.py) |
13+
14+
## Section 2.2 : Multiple Shooting Methods in Time
15+
16+
| 📜 Description | 🧮 Matlab | 🐍 Python |
17+
| :--- | :--- | :--- |
18+
| Multiple Shooting method | [MultipleShooting.m](./sec2.2/MultipleShooting.m), [CForwardEuler.m](./sec2.2/CForwardEuler.m) | [MultipleShooting.py](./sec2.2/MultipleShooting.py) |
19+
| example for Multiple Shooting on the Lorenz system | [exampleMultipleShootingLorenz.m](./sec2.2/exampleMultipleShootingLorenz.m) | [exampleMultipleShootingLorenz.py](./sec2.2/exampleMultipleShootingLorenz.py) |
20+
| scripts to generate the figures | [exampleMultipleShootingLorenzFig.m](./sec2.2/exampleMultipleShootingLorenzFig.m), [exampleMultipleShootingLorenzErrorFig.m](./sec2.2/exampleMultipleShootingLorenzErrorFig.m) | [figures.py](./sec2.2/figures.py) |
21+
22+
## Section 2.3 : The Parareal Algorithm
23+
24+
| 📜 Description | 🧮 Matlab | 🐍 Python |
25+
| :--- | :--- | :--- |
26+
| implementation of Parareal | [Parareal.m](./sec2.3/Parareal.m) | [Parareal.py](./sec2.3/Parareal.py) |
27+
| Parareal on the Lorenz equations | [examplePararealLorenz.m](./sec2.3/examplePararealLorenz.m) | [examplePararealLorenz.py](./sec2.3/examplePararealLorenz.py) |
28+
| error for Parareal on Lorenz equations | [examplePararealLorenzFig.m](./sec2.3/examplePararealLorenzFig.m) | [examplePararealLorenzFig.py](./sec2.3/examplePararealLorenzFig.py) |
29+
| Backward Euler solver for the Dahlquist equation | [DahlquistBE.m](./sec2.3/DahlquistBE.m), [SDahlquistBE.m](./sec2.3/SDahlquistBE.m) | [DahlquistBE.py](./sec2.3/DahlquistBE.py) |
30+
| Parareal on the Dahlquist equation | [examplePararealDahlquist.m](./sec2.3/examplePararealDahlquist.m) | [examplePararealDahlquist.py](./sec2.3/examplePararealDahlquist.py) |
31+
| Backward Euler solver for the Heat equation | [HeatEquationBE.m](./sec2.3/HeatEquationBE.m), [SHeatEquationBE.m](./sec2.3/SHeatEquationBE.m) | [HeatEquationBE.py](./sec2.3/HeatEquationBE.py) |
32+
| Parareal on the Heat equation | [examplePararealHeat.m](./sec2.3/examplePararealHeat.m) | [examplePararealHeat.py](./sec2.3/examplePararealHeat.py) |
33+
34+
[📖 Table of Content](../README.md)

‎chap2/sec2.1/FindClosestInterval.py

+27
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,27 @@
1+
def findClosestInterval(x,y):
2+
"""
3+
Finds the closest interval in a vector to a number.
4+
5+
Parameters
6+
----------
7+
x : float
8+
The value for which to find the closest interval.
9+
y : list or numpy.ndarray
10+
A sorted vector defining the intervals.
11+
12+
Returns
13+
-------
14+
int
15+
The lower index `m` of the interval in `y` that contains `x`.
16+
If `x` is outside the range of `y`, returns the closest interval.
17+
"""
18+
n = len(y)
19+
i = 0
20+
while i < n and x > y[i]:
21+
i += 1
22+
if i == 0:
23+
return 0
24+
elif i == n:
25+
return n-1
26+
else:
27+
return i-1

‎chap2/sec2.1/ForwardEuler.py

+5-3
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,8 @@ def forwardEuler(f, tspan, u0, N):
77
Parameters
88
----------
99
f : callable
10-
Function defining the system of ODEs, f(t, u).
10+
Function defining the system of ODEs, f(t, u),
11+
has to return a np.ndarray
1112
tspan : tuple
1213
Time interval as (t0, t_end).
1314
u0 : array-like
@@ -24,10 +25,11 @@ def forwardEuler(f, tspan, u0, N):
2425
"""
2526
dt = (tspan[1] - tspan[0]) / N # Time step size
2627
t = np.linspace(tspan[0], tspan[1], N + 1) # Time points
28+
u0 = np.asarray(u0)
2729
u = np.zeros((N + 1, len(u0))) # Solution array
2830
u[0, :] = u0 # Set initial condition
2931

30-
for n in range(N):
31-
u[n + 1, :] = u[n, :] + dt * f(t[n], u[n, :]) # Forward Euler step
32+
for n in range(N): # Forward Euler steps
33+
u[n + 1, :] = u[n, :] + dt * f(t[n], u[n, :])
3234

3335
return t, u

‎chap2/sec2.1/Nievergelt.py

+88
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,88 @@
1+
import numpy as np
2+
3+
def nievergelt(solver, T, u0, N, uPred, Mn, width):
4+
"""
5+
Implementation of Nievergelt's method for scalar ODEs.
6+
7+
Parameters
8+
----------
9+
solver : callable
10+
Function of the form solver(t0, t1, u0) that solves the ODE numerically
11+
between t0 and t1 with u0 as the initial value, and returns the solution
12+
over the interval [t0, t1].
13+
T : float
14+
The total time interval [0, T].
15+
u0 : float
16+
The initial value of the ODE.
17+
N : int
18+
The number of time sub-intervals dividing (0, T].
19+
uPred : ndarray
20+
Predicted solution, typically a cheap approximation.
21+
Mn : int
22+
The number of initial values for each sub-interval.
23+
width : float
24+
The width of the interval, centered around `U^0_n`, on which
25+
the initial values are uniformly distributed.
26+
27+
Returns
28+
-------
29+
U : ndarray
30+
The Nievergelt solution at the sub-interval endpoints.
31+
uTraj : ndarray
32+
The trajectories computed for each sub-interval.
33+
"""
34+
dT = T / N
35+
TT = np.linspace(0, T, N + 1)
36+
uTraj = np.zeros((N, Mn, len(solver(TT[0], TT[1], u0))))
37+
38+
# Compute first trajectory
39+
uTraj[0, 0, :] = solver(TT[0], TT[1], u0)
40+
41+
# Compute remaining trajectories
42+
for n in range(1, N):
43+
uStart = uPred[n] + np.linspace(-width / 2, width / 2, Mn)
44+
for m in range(Mn):
45+
uTraj[n, m, :] = solver(TT[n], TT[n + 1], uStart[m])
46+
47+
# Compute Nievergelt solution
48+
U = np.zeros(N + 1)
49+
U[0] = u0
50+
U[1] = uTraj[0, 0, -1]
51+
52+
for n in range(1, N):
53+
m = findClosestInterval(U[n], uTraj[n, :, 0])
54+
uS1 = uTraj[n, m, 0]
55+
uS2 = uTraj[n, m + 1, 0]
56+
p = (U[n] - uS2) / (uS1 - uS2)
57+
U[n + 1] = p * uTraj[n, m, -1] + (1 - p) * uTraj[n, m + 1, -1]
58+
59+
return U, uTraj
60+
61+
62+
def findClosestInterval(x,y):
63+
"""
64+
Finds the closest interval in a vector to a number.
65+
66+
Parameters
67+
----------
68+
x : float
69+
The value for which to find the closest interval.
70+
y : list or numpy.ndarray
71+
A sorted vector defining the intervals.
72+
73+
Returns
74+
-------
75+
int
76+
The lower index `m` of the interval in `y` that contains `x`.
77+
If `x` is outside the range of `y`, returns the closest interval.
78+
"""
79+
n = len(y)
80+
i = 0
81+
while i < n and x > y[i]:
82+
i += 1
83+
if i == 0:
84+
return 0
85+
elif i == n:
86+
return n-1
87+
else:
88+
return i-1

‎chap2/sec2.1/README.md

-10
This file was deleted.

‎chap2/sec2.1/exampleNievergelt.py

+12
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
import numpy as np
2+
from Nievergelt import nievergelt
3+
from ForwardEuler import forwardEuler
4+
5+
f = lambda t,u: np.cos(t)*u # RHS of the ODE problem
6+
u0 = 1; T = 2*np.pi # initial value, final time
7+
N = 10 # number of subintervals
8+
nSteps = 100 # fine steps per subinterval
9+
tPred, uPred = forwardEuler(f, [0, T], u0, N) # approximate prediction
10+
Mn = 2; width = 0.75 # algorithm parameters
11+
solver = lambda t0, t1, u0: forwardEuler(f, [t0, t1], u0, nSteps)[-1].ravel()
12+
U, uTraj = nievergelt(solver, T, u0, N, uPred, Mn, width)

‎chap2/sec2.1/exampleNievergeltLinearFig.m

+1-1
Original file line numberDiff line numberDiff line change
@@ -28,4 +28,4 @@
2828
end
2929
end
3030
legend(Location='northeast'), set(gca,'fontsize', 12);
31-
saveas(gcf(), '../Figures/NievergeltExampleLinear', 'epsc')
31+
saveas(gcf(), 'NievergeltExampleLinear', 'epsc')

‎chap2/sec2.1/exampleNievergeltLinearFig.py

+50
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
import numpy as np
2+
import matplotlib.pyplot as plt
3+
from Nievergelt import nievergelt
4+
from ForwardEuler import forwardEuler
5+
6+
# RHS of the ODE problem
7+
f = lambda t, u: np.cos(t) * u
8+
u0 = 1 # Initial solution
9+
T = 2 * np.pi # Final time
10+
N = 10 # Number of subintervals
11+
nSteps = 100 # Fine steps per subinterval
12+
13+
# Approximate prediction using Forward Euler
14+
tPred, uPred = forwardEuler(f, [0, T], u0, N)
15+
16+
# Nievergelt's method parameters
17+
Mn = 2
18+
width = 0.75
19+
20+
# Solver function for Nievergelt's method
21+
solver = lambda t0, t1, u0: forwardEuler(f, [t0, t1], u0, nSteps)[-1].ravel()
22+
23+
# Nievergelt's method
24+
U, uTraj = nievergelt(solver, T, u0, N, uPred, Mn, width)
25+
26+
# Fine solution for comparison
27+
tFine, uFine = forwardEuler(f, [0, T], u0, N * nSteps)
28+
29+
# Plotting
30+
plt.plot(tFine, uFine, label='Fine', zorder=3)
31+
plt.xlabel('Time')
32+
plt.ylabel('Solution')
33+
plt.plot(tPred, uPred, 'o', label='Prediction', zorder=2)
34+
plt.plot(tPred, U, '^', label='Nievergelt', zorder=2)
35+
36+
# Plot trajectories
37+
TT = np.linspace(0, T, N + 1)
38+
x = uTraj[0, 0, :]
39+
plt.plot(np.linspace(TT[0], TT[1], nSteps + 1), x, 'k--', label='Trajectories', zorder=1)
40+
41+
for n in range(1, N):
42+
for m in range(Mn):
43+
x = uTraj[n, m, :]
44+
plt.plot(np.linspace(TT[n], TT[n + 1], nSteps + 1), x, 'k--', zorder=1, alpha=0.7)
45+
46+
# Add legend and finalize plot
47+
plt.legend(loc='upper right')
48+
plt.gca().tick_params(labelsize=12)
49+
plt.savefig('NievergeltExampleLinear.eps', format='eps')
50+
plt.show()

‎chap2/sec2.1/exampleNievergeltNonLinearFig.m

+2-2
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@
2929
end
3030
end
3131
legend(Location='northeast'), set(gca,'fontsize', 12);
32-
saveas(figure(1), '../Figures/NievergeltExampleNonLinear_1', 'epsc')
32+
saveas(figure(1), 'NievergeltExampleNonLinear_1', 'epsc')
3333

3434
Mn=4; width=0.4; % Nievergelt's method parameters
3535
solver=@(t0,t1,u0) UForwardEuler(f,t0,t1,u0,nSteps);
@@ -53,4 +53,4 @@
5353
end
5454
end
5555
legend(Location='northeast'), set(gca,'fontsize', 12);
56-
saveas(figure(2), '../Figures/NievergeltExampleNonLinear_2', 'epsc')
56+
saveas(figure(2), 'NievergeltExampleNonLinear_2', 'epsc')

‎chap2/sec2.1/exampleNievergeltNonLinearFig.py

+76
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,76 @@
1+
import numpy as np
2+
import matplotlib.pyplot as plt
3+
from Nievergelt import nievergelt
4+
from ForwardEuler import forwardEuler
5+
6+
# Define the RHS of the ODE problem
7+
f = lambda t, u: np.cos(t) * np.exp(-u)
8+
u0 = np.log(2) # Initial solution
9+
T = 2 * np.pi # Final time
10+
N = 10 # Number of subintervals
11+
nSteps = 100 # Fine steps per subinterval
12+
13+
# Approximate prediction using Forward Euler
14+
tPred, uPred = forwardEuler(f, [0, T], u0, N)
15+
16+
# First set of Nievergelt's method parameters
17+
Mn = 2
18+
width = 0.2
19+
solver = lambda t0, t1, u0: forwardEuler(f, [t0, t1], u0, nSteps)[-1].ravel()
20+
U, uTraj = nievergelt(solver, T, u0, N, uPred, Mn, width)
21+
22+
# Fine solution for comparison
23+
tFine, uFine = forwardEuler(f, [0, T], u0, N * nSteps)
24+
25+
# Plot the first set of results
26+
plt.figure(1)
27+
plt.plot(tFine, uFine, label='Fine', zorder=3)
28+
plt.xlabel('Time')
29+
plt.ylabel('Solution')
30+
plt.plot(tPred, uPred, 'o', label='Prediction', zorder=2)
31+
plt.plot(tPred, U, '^', label='Nievergelt', zorder=2)
32+
33+
# Plot trajectories
34+
TT = np.linspace(0, T, N + 1)
35+
x = uTraj[0, 0, :]
36+
plt.plot(np.linspace(TT[0], TT[1], nSteps + 1), x, 'k--', label='Trajectories', zorder=1)
37+
38+
for n in range(1, N):
39+
for m in range(Mn):
40+
x = uTraj[n, m, :]
41+
plt.plot(np.linspace(TT[n], TT[n + 1], nSteps + 1), x, 'k--', zorder=1, alpha=0.7)
42+
43+
plt.legend(loc='upper right')
44+
plt.gca().tick_params(labelsize=12)
45+
plt.savefig('NievergeltExampleNonLinear_1.eps', format='eps')
46+
plt.show()
47+
48+
# Second set of Nievergelt's method parameters
49+
Mn = 4
50+
width = 0.4
51+
U, uTraj = nievergelt(solver, T, u0, N, uPred, Mn, width)
52+
53+
# Fine solution for comparison
54+
tFine, uFine = forwardEuler(f, [0, T], u0, N * nSteps)
55+
56+
# Plot the second set of results
57+
plt.figure(2)
58+
plt.plot(tFine, uFine, label='Fine', zorder=3)
59+
plt.xlabel('Time')
60+
plt.ylabel('Solution')
61+
plt.plot(tPred, uPred, 'o', label='Prediction', zorder=2)
62+
plt.plot(tPred, U, '^', label='Nievergelt', zorder=2)
63+
64+
# Plot trajectories
65+
x = uTraj[0, 0, :]
66+
plt.plot(np.linspace(TT[0], TT[1], nSteps + 1), x, 'k--', label='Trajectories', zorder=1)
67+
68+
for n in range(1, N):
69+
for m in range(Mn):
70+
x = uTraj[n, m, :]
71+
plt.plot(np.linspace(TT[n], TT[n + 1], nSteps + 1), x, 'k--', zorder=1, alpha=0.7)
72+
73+
plt.legend(loc='upper right')
74+
plt.gca().tick_params(labelsize=12)
75+
plt.savefig('NievergeltExampleNonLinear_2.eps', format='eps')
76+
plt.show()

0 commit comments

Comments
 (0)