2 Commits

Author SHA1 Message Date
327af48e73 Remove comment for problem 3 2023-09-01 16:22:54 +02:00
450e1859b1 Solve prblem 3 and add it to latex file 2023-09-01 16:21:27 +02:00
20 changed files with 105 additions and 464 deletions

8
.gitignore vendored
View File

@@ -36,11 +36,3 @@
*.log
*.out
*.bib
*.synctex.gz
*.bbl
# C++ specifics
src/*
!src/Makefile
!src/*.cpp
!src/*.py

0
include/problem2.hpp Normal file
View File

Binary file not shown.

View File

@@ -1,9 +1,7 @@
\documentclass[english,notitlepage]{revtex4-1} % defines the basic parameters of the document
%For preview: skriv i terminal: latexmk -pdf -pvc filnavn
% Silence warning of revtex4-1
\usepackage{silence}
\WarningFilter{revtex4-1}{Repair the float}
% if you want a single-column, remove reprint
@@ -15,7 +13,7 @@
%% I recommend downloading TeXMaker, because it includes a large library of the most common packages.
\usepackage{physics,amssymb} % mathematical symbols (physics imports amsmath)
\usepackage{amsmath}
\include{amsmath}
\usepackage{graphicx} % include graphics such as plots
\usepackage{xcolor} % set colors
\usepackage{hyperref} % automagic cross-referencing (this is GODLIKE)
@@ -74,9 +72,6 @@
%%
%% Don't ask me why, I don't know.
% custom stuff
\graphicspath{{./images/}}
\begin{document}
\title{Project 1} % self-explanatory
@@ -89,22 +84,116 @@
\textit{https://github.uio.no/FYS3150-G2-2023/Project-1}
\input{problems/problem1}
\section*{Problem 1}
\input{problems/problem2}
\begin{align*}
u(x) &= \int \int \frac{d^2 u}{dx^2} dx^2\\
&= \int \int -100 e^{-10x} dx^2 \\
&= \int \frac{-100 e^{-10x}}{-10} + c_1 dx \\
&= \int 10 e^{-10x} + c_1 dx \\
&= \frac{10 e^{-10x}}{-10} + c_1 x + c_2 \\
&= -e^{-10x} + c_1 x + c_2
\end{align*}
\input{problems/problem3}
Using the boundary conditions, we can find $c_1$ and $c_2$ as shown below:
\input{problems/problem4}
\begin{align*}
u(0) &= 0 \\
-e^{-10 \cdot 0} + c_1 \cdot 0 + c_2 &= 0 \\
-1 + c_2 &= 0 \\
c_2 &= 1
\end{align*}
\input{problems/problem5}
\begin{align*}
u(1) &= 0 \\
-e^{-10 \cdot 1} + c_1 \cdot 1 + c_2 &= 0 \\
-e^{-10} + c_1 + c_2 &= 0 \\
c_1 &= e^{-10} - c_2\\
c_1 &= e^{-10} - 1\\
\end{align*}
\input{problems/problem6}
Using the values that we found for $c_1$ and $c_2$, we get
\input{problems/problem7}
\begin{align*}
u(x) &= -e^{-10x} + (e^{-10} - 1) x + 1 \\
&= 1 - (1 - e^{-10}) - e^{-10x}
\end{align*}
\input{problems/problem8}
\section*{Problem 2}
\input{problems/problem9}
% Write which .cpp/.hpp/.py (using a link?) files are relevant for this and show the plot generated.
\section*{Problem 3}
To derive the discretized version of the Poisson equation, we first need
the taylor expansion for $u(x)$ around $x + h$ and $x - h$.
\begin{align*}
u(x+h) &= u(x) + u'(x) h + \frac{1}{2} u''(x) h^2 + \frac{1}{6} u'''(x) h^3 + \mathcal{O}(h^4)
\end{align*}
\begin{align*}
u(x-h) &= u(x) - u'(x) h + \frac{1}{2} u''(x) h^2 - \frac{1}{6} u'''(x) h^3 + \mathcal{O}(h^4)
\end{align*}
If we add the equations above, we get this new equation:
\begin{align*}
u(x+h) + u(x-h) &= 2 u(x) + u''(x) h^2 + \mathcal{O}(h^4) \\
u(x+h) - 2 u(x) + u(x-h) + \mathcal{O}(h^4) &= u''(x) h^2 \\
u''(x) &= \frac{u(x+h) - 2 u(x) + u(x-h)}{h^2} + \mathcal{O}(h^2) \\
u_i''(x) &= \frac{u_{i+1} - 2 u_i + u_{i-1}}{h^2} + \mathcal{O}(h^2) \\
\end{align*}
We can then replace $\frac{d^2u}{dx^2}$ with the RHS (right-hand side) of the equation:
\begin{align*}
- \frac{d^2u}{dx^2} &= 100 e^{-10x} \\
\frac{ - u_{i+1} + 2 u_i - u_{i-1}}{h^2} + \mathcal{O}(h^2) &= 100 e^{-10x} \\
\end{align*}
And lastly, we leave out $\mathcal{O}(h^2)$ and change $u_i$ to $v_i$ to
differentiate between the exact solution and the approximate solution,
and get the discretized version of the equation:
\begin{align*}
align* \frac{ - u_{i+1} + 2 u_i - u_{i-1}}{h^2} &= 100 e^{-10x} \\
\end{align*}
\section*{Problem 4}
% Show that each iteration of the discretized version naturally creates a matrix equation.
\section*{Problem 5}
\subsection*{a)}
\subsection*{b)}
\section*{Problem 6}
\subsection*{a)}
% Use Gaussian elimination, and then use backwards substitution to solve the equation
\subsection*{b)}
% Figure it out
\section*{Problem 7}
% Link to relevant files on gh and possibly add some comments
\section*{Problem 8}
%link to relvant files and show plots
\section*{Problem 9}
% Show the algorithm, then calculate FLOPs, then link to relevant files
\section*{Problem 10}
% Time and show result, and link to relevant files
\end{document}

Binary file not shown.

View File

@@ -1,44 +0,0 @@
\section*{Problem 1}
First, we rearrange the equation.
\begin{align*}
- \frac{d^2u}{dx^2} &= 100 e^{-10x} \\
\frac{d^2u}{dx^2} &= -100 e^{-10x} \\
\end{align*}
Now we find $u(x)$.
% Do the double integral
\begin{align*}
u(x) &= \int \int \frac{d^2 u}{dx^2} dx^2 \\
&= \int \int -100 e^{-10x} dx^2 \\
&= \int \frac{-100 e^{-10x}}{-10} + c_1 dx \\
&= \int 10 e^{-10x} + c_1 dx \\
&= \frac{10 e^{-10x}}{-10} + c_1 x + c_2 \\
&= -e^{-10x} + c_1 x + c_2
\end{align*}
Using the boundary conditions, we can find $c_1$ and $c_2$
\begin{align*}
u(0) &= 0 \\
-e^{-10 \cdot 0} + c_1 \cdot 0 + c_2 &= 0 \\
-1 + c_2 &= 0 \\
c_2 &= 1
\end{align*}
\begin{align*}
u(1) &= 0 \\
-e^{-10 \cdot 1} + c_1 \cdot 1 + c_2 &= 0 \\
-e^{-10} + c_1 + c_2 &= 0 \\
c_1 &= e^{-10} - c_2\\
c_1 &= e^{-10} - 1\\
\end{align*}
Using the values that we found for $c_1$ and $c_2$, we get
\begin{align*}
u(x) &= -e^{-10x} + (e^{-10} - 1) x + 1 \\
&= 1 - (1 - e^{-10}) - e^{-10x}
\end{align*}

View File

@@ -1,3 +0,0 @@
\section*{Problem 10}
% Time and show result, and link to relevant files

View File

@@ -1,11 +0,0 @@
\section*{Problem 2}
The code for generating the points and plotting them can be found under.
Point generator code: https://github.uio.no/FYS3150-G2-203/Project-1/blob/main/src/analyticPlot.cpp
Plotting code: https://github.uio.no/FYS3150-G2-2023/Project-1/blob/main/src/analyticPlot.py
Here is the plot of the analytical solution for $u(x)$.
\includegraphics[scale=.5]{analytical_solution.pdf}

View File

@@ -1,37 +0,0 @@
\section*{Problem 3}
To derive the discretized version of the Poisson equation, we first need
the Taylor expansion for $u(x)$ around $x$ for $x + h$ and $x - h$.
\begin{align*}
u(x+h) &= u(x) + u'(x) h + \frac{1}{2} u''(x) h^2 + \frac{1}{6} u'''(x) h^3 + \mathcal{O}(h^4)
\end{align*}
\begin{align*}
u(x-h) &= u(x) - u'(x) h + \frac{1}{2} u''(x) h^2 - \frac{1}{6} u'''(x) h^3 + \mathcal{O}(h^4)
\end{align*}
If we add the equations above, we get this new equation:
\begin{align*}
u(x+h) + u(x-h) &= 2 u(x) + u''(x) h^2 + \mathcal{O}(h^4) \\
u(x+h) - 2 u(x) + u(x-h) + \mathcal{O}(h^4) &= u''(x) h^2 \\
u''(x) &= \frac{u(x+h) - 2 u(x) + u(x-h)}{h^2} + \mathcal{O}(h^2) \\
u_i''(x) &= \frac{u_{i+1} - 2 u_i + u_{i-1}}{h^2} + \mathcal{O}(h^2) \\
\end{align*}
We can then replace $\frac{d^2u}{dx^2}$ with the RHS (right-hand side) of the equation:
\begin{align*}
- \frac{d^2u}{dx^2} &= f(x) \\
\frac{ - u_{i+1} + 2 u_i - u_{i-1}}{h^2} + \mathcal{O}(h^2) &= f_i \\
\end{align*}
And lastly, we leave out $\mathcal{O}(h^2)$ and change $u_i$ to $v_i$ to
differentiate between the exact solution and the approximate solution,
and get the discretized version of the equation:
\begin{align*}
\frac{ - v_{i+1} + 2 v_i - v_{i-1}}{h^2} &= 100 e^{-10x_i} \\
\end{align*}

View File

@@ -1,44 +0,0 @@
\section*{Problem 4}
% Show that each iteration of the discretized version naturally creates a matrix equation.
The value of $u(x_{0})$ and $u(x_{1})$ is known, using the discretized equation we can approximate the value of $f(x_{i}) = f_{i}$. This will result in a set of equations
\begin{align*}
- v_{0} + 2 v_{1} - v_{2} &= h^{2} \cdot f_{1} \\
- v_{1} + 2 v_{2} - v_{3} &= h^{2} \cdot f_{2} \\
\vdots & \\
- v_{m-2} + 2 v_{m-1} - v_{m} &= h^{2} \cdot f_{m-1} \\
\end{align*}
Rearranging the first and last equation, moving terms of known boundary values to the RHS
\begin{align*}
2 v_{1} - v_{2} &= h^{2} \cdot f_{1} + v_{0} \\
- v_{1} + 2 v_{2} - v_{3} &= h^{2} \cdot f_{2} \\
\vdots & \\
- v_{m-2} + 2 v_{m-1} &= h^{2} \cdot f_{m-1} + v_{m} \\
\end{align*}
We now have a number of linear eqations, corresponding to the number of unknown values, which can be represented as an augmented matrix
\begin{align*}
\left[
\begin{matrix}
2v_{1} & -v_{2} & 0 & \dots & 0 \\
-v_{1} & 2v_{2} & -v_{3} & 0 & \\
0 & -v_{2} & 2v_{3} & -v_{4} & \\
\vdots & & & \ddots & \vdots \\
0 & & & -v_{m-2} & 2v_{m-1} \\
\end{matrix}
\left|
\,
\begin{matrix}
g_{1} \\
g_{2} \\
g_{2} \\
\vdots \\
g_{m-1} \\
\end{matrix}
\right.
\right]
\end{align*}
where $g_{i} = h^{2} f_{i}$. An augmented matrix can be represented as $\boldsymbol{A} \vec{x} = \vec{b}$. In this case $\boldsymbol{A}$ is the coefficient matrix with a tridiagonal signature $(-1, 2, -1)$ and dimension $n \cross n$, where $n=m-2$.

View File

@@ -1,6 +0,0 @@
\section*{Problem 5}
\subsection*{a)}
\subsection*{b)}

View File

@@ -1,9 +0,0 @@
\section*{Problem 6}
\subsection*{a)}
% Use Gaussian elimination, and then use backwards substitution to solve the equation
\subsection*{b)}
% Figure it out

View File

@@ -1,3 +0,0 @@
\section*{Problem 7}
% Link to relevant files on gh and possibly add some comments

View File

@@ -1,3 +0,0 @@
\section*{Problem 8}
%link to relvant files and show plots

View File

@@ -1,3 +0,0 @@
\section*{Problem 9}
% Show the algorithm, then calculate FLOPs, then link to relevant files

View File

@@ -1,17 +0,0 @@
CC=g++
.PHONY: clean
all: simpleFile analyticPlot
simpleFile: simpleFile.o
$(CC) -o $@ $^
analyticPlot: analyticPlot.o
$(CC) -o $@ $^
%.o: %.cpp
$(CC) -c $< -o $@
clean:
rm *.o

View File

@@ -1,55 +0,0 @@
#include <iostream>
#include <cmath>
#include <vector>
#include <string>
#include <numeric>
#include <fstream>
#include <iomanip>
#define RANGE 1000
#define FILENAME "analytical_solution.txt"
double u(double x);
void generate_range(std::vector<double> &vec, double start, double stop, int n);
void write_analytical_solution(std::string filename, int n);
int main() {
write_analytical_solution(FILENAME, RANGE);
return 0;
};
double u(double x) {
return 1 - (1 - exp(-10))*x - exp(-10*x);
};
void generate_range(std::vector<double> &vec, double start, double stop, int n) {
double step = (stop - start) / n;
for (int i = 0; i <= vec.size(); i++) {
vec[i] = i * step;
}
}
void write_analytical_solution(std::string filename, int n) {
std::vector<double> x(n), y(n);
generate_range(x, 0.0, 1.0, n);
// Set up output file and strem
std::ofstream outfile;
outfile.open(filename);
// Parameters for formatting
int width = 12;
int prec = 4;
// Calculate u(x) and write to file
for (int i = 0; i <= x.size(); i++) {
y[i] = u(x[i]);
outfile << std::setw(width) << std::setprecision(prec) << std::scientific << x[i]
<< std::setw(width) << std::setprecision(prec) << std::scientific << y[i]
<< std::endl;
}
outfile.close();
}

View File

@@ -1,19 +0,0 @@
import numpy as np
import matplotlib.pyplot as plt
def main():
FILENAME = "analytical_solution.pdf"
x = []
v = []
with open('analytical_solution.txt') as f:
for line in f:
a, b = line.strip().split()
x.append(float(a))
v.append(float(b))
plt.plot(x, v)
plt.savefig(FILENAME)
if __name__ == "__main__":
main()

View File

@@ -1,24 +0,0 @@
#include "GeneralAlgorithm.hpp"
#include <armadillo>
#include <iostream>
double f(double x) {
return 100. * std::exp(-10.*x);
}
double a_sol(double x) {
return 1. - (1. - std::exp(-10)) * x - std::exp(-10*x);
}
int main() {
arma::mat A = arma::eye(3,3);
GeneralAlgorithm ga(3, &A, f, a_sol, 0., 1.);
ga.solve();
std::cout << "Time: " << ga.time(5) << std::endl;
ga.error();
return 0;
}

View File

@@ -1,162 +0,0 @@
#include <armadillo>
#include <cmath>
#include <ctime>
#include <fstream>
#include <iomanip>
#include <ios>
#include <string>
#define TIMING_ITERATIONS 5
arma::vec* general_algorithm(
arma::vec* sub_diag,
arma::vec* main_diag,
arma::vec* sup_diag,
arma::vec* g_vec
)
{
int n = g_vec->n_elem;
double d;
for (int i = 1; i < n; i++) {
d = (*sub_diag)(i-1) / (*main_diag)(i-1);
(*main_diag)(i) -= d*(*sup_diag)(i-1);
(*g_vec)(i) -= d*(*g_vec)(i-1);
}
(*g_vec)(n-1) /= (*main_diag)(n-1);
for (int i = n-2; i >= 0; i--) {
(*g_vec)(i) = ((*g_vec)(i) - (*sup_diag)(i) * (*g_vec)(i+1)) / (*main_diag)(i);
}
return g_vec;
}
arma::vec* special_algorithm(
double sub_sig,
double main_sig,
double sup_sig,
arma::vec* g_vec
)
{
int n = g_vec->n_elem;
arma::vec diag = arma::vec(n);
for (int i = 1; i < n; i++) {
// Calculate values for main diagonal based on indices
diag(i-1) = (double)(i+1) / i;
(*g_vec)(i) += (*g_vec)(i-1) / diag(i-1);
}
// The last element in main diagonal has value (i+1)/i = (n+1)/n
(*g_vec)(n-1) /= (double)(n+1) / (n);
for (int i = n-2; i >= 0; i--) {
(*g_vec)(i) = ((*g_vec)(i) + (*g_vec)(i+1))/ diag(i);
}
return g_vec;
}
void error(
std::string filename,
arma::vec* x_vec,
arma::vec* v_vec,
arma::vec* a_vec
)
{
std::ofstream ofile;
ofile.open(filename);
if (!ofile.is_open()) {
exit(1);
}
for (int i=0; i < a_vec->n_elem; i++) {
double sub = (*a_vec)(i) - (*v_vec)(i);
ofile << std::setprecision(8) << std::scientific << (*x_vec)(i)
<< std::setprecision(8) << std::scientific << std::log10(std::abs(sub))
<< std::setprecision(8) << std::scientific << std::log10(std::abs(sub/(*a_vec)(i)))
<< std::endl;
}
ofile.close();
}
double f(double x) {
return 100*std::exp(-10*x);
}
void build_array(
int n_steps,
arma::vec* sub_diag,
arma::vec* main_diag,
arma::vec* sup_diag,
arma::vec* g_vec
)
{
sub_diag->resize(n_steps-2);
main_diag->resize(n_steps-1);
sup_diag->resize(n_steps-2);
sub_diag->fill(-1);
main_diag->fill(2);
sup_diag->fill(-1);
g_vec->resize(n_steps-1);
double step_size = 1./ (double) n_steps;
for (int i=0; i < n_steps-1; i++) {
(*g_vec)(i) = f((i+1)*step_size);
}
}
void timing() {
arma::vec sub_diag, main_diag, sup_diag, g_vec;
int n_steps;
std::ofstream ofile;
ofile.open("timing.txt");
// Timing
for (int i=1; i <= 8; i++) {
n_steps = std::pow(10, i);
clock_t g_1, g_2, s_1, s_2;
double g_res = 0, s_res = 0;
for (int j=0; j < TIMING_ITERATIONS; j++) {
build_array(n_steps, &sub_diag, &main_diag, &sup_diag, &g_vec);
g_1 = clock();
general_algorithm(&sub_diag, &main_diag, &sup_diag, &g_vec);
g_2 = clock();
g_res += (double) (g_2 - g_1) / CLOCKS_PER_SEC;
build_array(n_steps, &sub_diag, &main_diag, &sup_diag, &g_vec);
s_1 = clock();
special_algorithm(-1., 2., -1., &g_vec);
s_2 = clock();
s_res += (double) (s_2 - s_1) / CLOCKS_PER_SEC;
}
ofile
<< n_steps << ","
<< g_res / (double) TIMING_ITERATIONS << ","
<< s_res / (double) TIMING_ITERATIONS << std::endl;
}
ofile.close();
}
int main()
{
timing();
}