lars/FYS4150
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fix issues with limited numerical precision and artifacts with big N.

Browse Source
This commit is contained in:
Lars Bogner
2025-09-04 08:56:55 +02:00
parent dc5fba552e
commit e62c42bb2c
5 changed files with 94 additions and 32 deletions
Download Patch File Download Diff File Expand all files Collapse all files
projects/project1/main.pdf
BIN
View File
Binary file not shown.
projects/project1/main.tex
+3 -2
View File
@@ -243,7 +243,7 @@ Plotting the absolute error for different values of $N$ reveals the said converg
\end{figure}
\subsection*{Subproblem (b)}
Plotting the relative error for different values of $N$ reveals the increase of relative error as $u$ approaches 0 in the vicinity of the boundaries. This behavior is prevalent in all solutions as the limit of small divisors is approached. This behavior is quite logical, but it also shows, that not all uncertainties can be correctly minimized with an increase in resolution. The results are shown in \autoref{fig:rel_error_plot}. To overcome the division by zero issue for the relative error, only the values of $u > \num{1e-10}$ are considered.
Plotting the relative error for different values of $N$ reveals a constant value for the relative error accross the entire domain of $x$. The relative error decreases by two orders of magnitude for every order of magnitude increase in $N$. The results are summarized in \autoref{fig:rel_error_plot}. To avoiud issues with divisions by zero, the points where $u(x) \leq \num{e-10}$ were omitted from the plot.
\begin{figure}
\centering
@@ -253,7 +253,7 @@ Plotting the relative error for different values of $N$ reveals the increase of
\end{figure}
\subsection*{Subproblem (c)}
The maximum relative errors are calculated directly using C++ as the output size scales linearly with $N$ and gets overwhelmingly large for bigger $N$. The results are output using the \texttt{std::cout} command and retrieved using GNU/Linux command line utilities. The code for calculating the maximum relative error $\max_i \frac{\delta v_i}{u_i}$ (see \eqref{eq:relerror}) is included in the C++ implementation. The results are summarized in \autoref{tab:rel_errors}. We see a steady decrease of the relative error up to $N = \num{100000}$, after which the relative error increases again. This will most likely be due to the artifacts seen in \autoref{fig:rel_error_plot} starting at $N = \num{10000}$. The source of these artifacts is not entirely clear.
The maximum relative errors are calculated directly using C++ as the output size scales linearly with $N$ and gets overwhelmingly large for bigger $N$. The results are output using the \texttt{std::cout} command and retrieved using GNU/Linux command line utilities. The code for calculating the maximum relative error $\max_i \frac{\delta v_i}{u_i}$ (see \eqref{eq:relerror}) is included in the C++ implementation. The results are summarized in \autoref{tab:rel_errors}. We see a steady decrease of the relative error up to $N = \num{100000}$, after which the relative error increases again. This will most likely be due to the limited numerical precision of floating point numbers in C++ and associated numerical errors.
\begin{table}[H]
\centering
@@ -332,3 +332,4 @@ To evaluate the performance of the different algorithms, the algorithms were exe
\bibliography{include/citations}
\end{document}
src/project1/poisson_solver.cpp
+4 -4
View File
@@ -43,11 +43,11 @@ int main(int argc, char* argv[]){
// Output the solution
string filename = "numerical_solution_" + to_string(N) + ".csv";
ofstream ofs(filename);
ofs << "x,u(x)" << endl;
ofs << "x,v(x),u(x)" << endl;
for (int i = 0; i <= N; i++) {
ofs << fixed << setprecision(decimal_places) << scientific << x[i] << "," << solution[i] << endl;
}
ofs.close();
ofs << fixed << setprecision(decimal_places) << scientific << x[i] << "," << solution[i] << "," << u[i] << endl;
}
ofs.close();
}
src/project1/python/poisson_plotter.py
+57 -26
View File
@@ -13,12 +13,43 @@ args = parser.parse_args()
FILENAME = args.filename
REFERENCE = args.reference
def get_rc_params():
colors = ["FF220C", "70D6FF", "8AAA79", "666370", "1C1F33"]
rcParams = plt.rcParams
# Use LaTeX for rendering
# Setup fonts
rcParams["text.usetex"] = True
rcParams["font.family"] = "serif"
rcParams["font.size"] = 12
rcParams["axes.labelsize"] = 12
rcParams["axes.titlesize"] = 12
rcParams["legend.fontsize"] = 10
rcParams["xtick.labelsize"] = 10
rcParams["ytick.labelsize"] = 10
# Figure size and resolution
rcParams["figure.figsize"] = (6, 4)
rcParams["figure.dpi"] = 300
# Use colors from the palette
rcParams["axes.prop_cycle"] = plt.cycler(color=[f"#{color}" for color in colors])
# Grid
rcParams["axes.grid"] = True
rcParams["grid.alpha"] = 0.5
rcParams["grid.linestyle"] = "--"
# Point ticks to the inside of the axes
rcParams["xtick.direction"] = "in"
rcParams["ytick.direction"] = "in"
rcParams["xtick.top"] = True
rcParams["ytick.right"] = True
return rcParams
def singular_plot():
plt.rcParams.update(get_rc_params())
df = pd.read_csv(FILENAME)
fig, ax = plt.subplots(figsize=(6,3))
ax.plot(df["x"], df["u(x)"])
ax.plot(df["x"], df[df.columns[1]])
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$u(x)$")
ax.set_ylabel(rf"${df.columns[1]}$")
fig.tight_layout()
@@ -27,22 +58,21 @@ def singular_plot():
fig.savefig(os.path.join(output_directory, output_file))
def reference_plot():
plt.rcParams.update(get_rc_params())
df = pd.read_csv(FILENAME).sort_values("x")
df_ref = pd.read_csv(REFERENCE).sort_values("x")
fig, (ax1, ax2) = plt.subplots(nrows=2, sharex=True, height_ratios=[3, 1], figsize=(4, 8))
ax1.plot(df["x"], df["u(x)"], label="Numerical Solution")
ax1.plot(df_ref["x"], df_ref["u(x)"], label="Reference Solution", linestyle="--")
ax1.plot(df["x"], df["v(x)"], label="Numerical Solution")
ax1.plot(df["x"], df["u(x)"], label="Reference Solution", linestyle="--")
ax1.set_xlabel(r"$x$")
ax1.set_ylabel(r"$u(x)$")
ax1.legend()
combined_df = pd.merge_asof(df, df_ref, on="x", suffixes=("_num", "_ref"))
ratio = (combined_df["u(x)_num"] - combined_df["u(x)_ref"]) / combined_df["u(x)_ref"]
ratio = (df["v(x)"] - df["u(x)"]) / df["u(x)"]
ax2.plot(combined_df["x"], ratio.abs())
ax2.plot(df["x"], ratio.abs())
ax2.set_xlabel(r"$x$")
ax2.set_ylabel(r"$|\frac{u(x)_{num} - u(x)_{ref}}{u(x)_{ref}}|$")
ax2.set_ylabel(r"$|\frac{v(x) - u(x)}{u(x)}|$")
ax2.axhline(y=1, color="gray", linestyle="--")
ax2.set_yscale("log")
@@ -53,30 +83,29 @@ def reference_plot():
fig.savefig(os.path.join(output_directory, output_file))
def multiple_reference_plot():
plt.rcParams.update(get_rc_params())
# Filename is wildcard -> Find files
import glob
files = glob.glob(FILENAME)
files.sort()
fig, (ax1, ax2) = plt.subplots(nrows=2, sharex=True, height_ratios=[3, 1], figsize=(6,5))
df_ref = pd.read_csv(REFERENCE).sort_values("x")
for file in files:
df = pd.read_csv(file).sort_values("x")
label = f"Num. (N = {os.path.basename(file).replace(".csv", "").split("_")[-1]})"
ax1.plot(df["x"], df["u(x)"], label=label)
ax1.plot(df["x"], df["v(x)"], label=label)
combined_df = pd.merge_asof(df, df_ref, on="x", suffixes=("_num", "_ref"))
ratio = (combined_df["u(x)_num"] - combined_df["u(x)_ref"]) / combined_df["u(x)_ref"]
ax2.plot(combined_df["x"], ratio.abs(), label=label)
ax1.plot(df_ref["x"], df_ref["u(x)"], label="Reference Solution", linestyle="--", color="gray")
ratio = (df["v(x)"] - df["u(x)"]) / df["u(x)"]
ax2.plot(df["x"], ratio.abs(), label=label)
ax1.plot(df["x"], df["u(x)"], label="Reference Solution", linestyle="--", color="gray")
ax1.set_xlabel(r"$x$")
ax1.set_ylabel(r"$u(x)$")
ax1.legend()
ax2.axhline(y=1, color="gray", linestyle="--")
ax2.set_xlabel(r"$x$")
ax2.set_ylabel(r"$|\frac{u(x)_{num} - u(x)_{ref}}{u(x)_{ref}}|$")
ax2.set_ylabel(r"$|\frac{v(x) - u(x)}{u(x)}|$")
ax2.set_yscale("log")
fig.tight_layout()
@@ -92,6 +121,7 @@ def error_plots():
def abs_error_plot():
plt.rcParams.update(get_rc_params())
import glob
files = glob.glob(FILENAME)
files.sort()
@@ -101,12 +131,12 @@ def abs_error_plot():
fig, ax = plt.subplots(figsize=(6,3))
for file in files:
df = pd.read_csv(file).sort_values("x")
combined_df = pd.merge_asof(df, df_ref, on="x", suffixes=("_num", "_ref"))
log_abs_err = np.log10((combined_df["u(x)_num"] - combined_df["u(x)_ref"]).abs())
ax.plot(combined_df["x"], log_abs_err, label=f"N={os.path.basename(file).replace('.csv','').split('_')[-1]}")
# combined_df = pd.merge_asof(df, df_ref, on="x", suffixes=("_num", "_ref"))
log_abs_err = np.log10((df["u(x)"] - df["v(x)"]).abs())
ax.plot(df["x"], log_abs_err, label=f"N={os.path.basename(file).replace('.csv','').split('_')[-1]}")
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$\log_{10}(|u(x)_{num} - u(x)_{ref}|)$")
ax.set_ylabel(r"$\log_{10}(|u(x) - v(x)|)$")
ax.legend()
fig.tight_layout()
@@ -117,6 +147,7 @@ def abs_error_plot():
def rel_error_plot():
plt.rcParams.update(get_rc_params())
import glob
files = glob.glob(FILENAME)
files.sort()
@@ -126,14 +157,14 @@ def rel_error_plot():
fig, ax = plt.subplots(figsize=(6,3))
for file in files:
df = pd.read_csv(file).sort_values("x")
combined_df = pd.merge_asof(df, df_ref, on="x", suffixes=("_num", "_ref"))
# combined_df = pd.merge_asof(df, df_ref, on="x", suffixes=("_num", "_ref"))
# Drop rows where u(x)_ref is zero
combined_df = combined_df[combined_df["u(x)_ref"].abs() > 1e-10]
log_rel_err = np.log10((combined_df["u(x)_num"] - combined_df["u(x)_ref"]).abs() / combined_df["u(x)_ref"].abs())
ax.plot(combined_df["x"], log_rel_err, label=f"N={os.path.basename(file).replace('.csv','').split('_')[-1]}")
df = df[df["u(x)"] > 1e-10]
log_rel_err = np.log10((df["u(x)"] - df["v(x)"]).abs() / df["u(x)"].abs())
ax.plot(df["x"], log_rel_err, label=f"N={os.path.basename(file).replace('.csv','').split('_')[-1]}")
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$\log_{10}(|\frac{u(x)_{num} - u(x)_{ref}}{u(x)_{ref}}|)$")
ax.set_ylabel(r"$\log_{10}(|\frac{u(x) - v(x)}{u(x)}|)$")
ax.legend()
fig.tight_layout()
src/project1/python/timing_plotter.py
+30
View File
@@ -4,10 +4,40 @@ import matplotlib.pyplot as plt
import argparse
import os
def get_rc_params():
colors = ["FF220C", "70D6FF", "8AAA79", "666370", "1C1F33"]
rcParams = plt.rcParams
# Use LaTeX for rendering
# Setup fonts
rcParams["text.usetex"] = True
rcParams["font.family"] = "serif"
rcParams["font.size"] = 12
rcParams["axes.labelsize"] = 12
rcParams["axes.titlesize"] = 12
rcParams["legend.fontsize"] = 10
rcParams["xtick.labelsize"] = 10
rcParams["ytick.labelsize"] = 10
# Figure size and resolution
rcParams["figure.figsize"] = (6, 4)
rcParams["figure.dpi"] = 300
# Use colors from the palette
rcParams["axes.prop_cycle"] = plt.cycler(color=[f"#{color}" for color in colors])
# Grid
rcParams["axes.grid"] = True
rcParams["grid.alpha"] = 0.5
rcParams["grid.linestyle"] = "--"
# Point ticks to the inside of the axes
rcParams["xtick.direction"] = "in"
rcParams["ytick.direction"] = "in"
rcParams["xtick.top"] = True
rcParams["ytick.right"] = True
return rcParams
parser = argparse.ArgumentParser(description="Plot timing results from CSV file.")
parser.add_argument("csv_file", help="Path to the CSV file containing timing results.")
args = parser.parse_args()
plt.rcParams.update(get_rc_params())
df = pd.read_csv(args.csv_file)
for column in df.columns[1:]:
df[column] = df[column] / df["N"] # Normalize by N