finished Ex3 and cleanup

This commit is contained in:
Daniel Kapla 2022-05-17 18:38:52 +02:00
parent c05ea06703
commit af315ca800
52 changed files with 434 additions and 652 deletions

1
.gitattributes vendored
View File

@ -1,2 +1,3 @@
*.png filter=lfs diff=lfs merge=lfs -text *.png filter=lfs diff=lfs merge=lfs -text
*.svg filter=lfs diff=lfs merge=lfs -text *.svg filter=lfs diff=lfs merge=lfs -text
*.mp4 filter=lfs diff=lfs merge=lfs -text

3
.gitignore vendored
View File

@ -10,3 +10,6 @@ plots/
# vscode workspace config # vscode workspace config
.vscode/ .vscode/
# LaTeX compile files

View File

@ -1,59 +0,0 @@
function [x,t,u] = Dirichlet_EA(L,N,T,K,c1,c2,f,u0,D)
% ---- Solution of heat equation ----
% u_t - u_xx = f in the interval [-L,L] (doesn't matter if we change it)
% Dirichlet BC
%
% -----------------------------------------------
% Sintassi:
% [x,t,u]=calore_template(L,N,T,K,c1,c2,fun,u0)
%
% Input:
% L Half of the width (-L,L)
% N number of intervals in (-L,L)
% T max time (0,T)
% K number of intervals in (0,T)
% c1 Dirichlet BC in x=-L
% c2 Dirichlet BC in x=L
% f force function
% u0 initial condition in t=0
%
% Output:
% x vector of the spatial nodes
% t vector of the time nodes
% u numeric solution of the problem
% discretisation step in time and space
h=L/N; %space
tau=T/K; %time
% Initialization of t
t=linspace(0,T,K+1)';
% initialization of x
x=linspace(0,L,N+1);
% Initialization of the matrix solution u
u=zeros(N+1,K+1);
% Intial conditions
u(:,1)=u0(x);
% BC
u(1,:)=c1(t);
u(end,:)=c2(t);
% Creation of the matrix A
% We obtained this matrix using the Central Discretization method
e=ones(N-1,1);
A=spdiags([-e,2*e,-e],[-1,0,1],N-1,N-1)/(h^2);
I=speye(N-1,N-1);
% Compute the solution
% What we are doing here is to compute the solution in the interval for
% each time k
for k=1:K
% Assembly of the force
F=f(x(2:end-1),t(k));
% Correction of the force using the BC
F(1)=F(1) + D*c1(t(k))/(h^2);
F(end)=F(end) + D*c2(t(k))/(h^2);
% Solution using Eulero forward
u(2:end-1,k+1) = ((I - tau*D*A)*u(2:end-1,k))' + tau*F;
end

Binary file not shown.

View File

@ -1,63 +0,0 @@
function [x,t,u] = Mixed_EA(L,N,T,K,c1,c2,f,u0,D)
% ---- Solution of heat equation ----
% u_t - u_xx = f in the interval [-L,L] (doesn't matter if we change it)
% Dirichlet BC
%
% -----------------------------------------------
% Sintassi:
% [x,t,u]=calore_template(L,N,T,K,c1,c2,fun,u0)
%
% Input:
% L Half of the width (-L,L)
% N number of intervals in (-L,L)
% T max time (0,T)
% K number of intervals in (0,T)
% c1 Dirichlet BC in x=-L
% c2 Dirichlet BC in x=L
% f force function
% u0 initial condition in t=0
%
% Output:
% x vector of the spatial nodes
% t vector of the time nodes
% u numeric solution of the problem
% discretisation step in time and space
h=L/N; %space
tau=T/K; %time
% Initialization of t
t=linspace(0,T,K+1)';
% initialization of x
x=linspace(0,L,N+1);
% Initialization of the matrix solution u
u=zeros(N+1,K+1);
% Intial conditions
u(:,1)=u0(x);
% BC Dirichlet
u(1,:)=c1(t);
% Creation of the matrix A
% We obtained this matrix using the Central Discretization method
e=ones(N-1,1);
A=spdiags([-e,2*e,-e],[-1,0,1],N-1,N-1)/(h^2);
%modify the matrix such that we can use the noimann condition to calculate
%the last node.
%In this case we want to use a second order decentralized approximation
A(end,end-1)=-2/(h^2);
I=speye(N-1,N-1);
% Compute the solution
% What we are doing here is to compute the solution in the interval for
% each time k
for k=1:K
% Assembly of the force
F=f(x(2:end-1),t(k));
% Correction of the force using the BC
F(1)=F(1) + D*c1(t(k))/(h^2);
F(end)=F(end) + 2*h*c2;
% Solution using Eulero forward
u(2:end-1,k+1) = ((I - D*tau*A)*u(2:end-1,k))' + tau*F;
u(end,k+1)=u(end-2,k+1) + c2*2*h*D;
end

View File

@ -1,62 +0,0 @@
function [x,t,u] = Mixed_EI(L,N,T,K,c1,c2,f,u0,D)
% ---- Risoluzione dell'equazione del calore ----
% u_t - u_xx = f nell'intervallo [-L,L]
% con condizioni al bordo di Dirichlet
% e condizioni iniziali.
% -----------------------------------------------
% Sintassi:
% [x,t,u]=calore_template(L,N,T,K,c1,c2,fun,u0)
%
% Input:
% L semiampiezza intervallo spaziale (-L,L)
% N numero di sottointervalli in (-L,L)
% T estremo finale intervallo temporale (0,T)
% K numero di sottointervalli in (0,T)
% c1 funzione che descrive la condizione di Dirichlet in x=-L
% c2 funzione che descrive la condizione di Dirichlet in x=L
% f funzione che descrive il termine noto dell'equazione
% u0 funzione che descrive la condizione iniziale in t=0
%
% Output:
% x vettore dei nodi spaziali
% t vettore dei nodi temporali
% u soluzione numerica
% del problema
% Calcolo passo di discretizzazione in spazio e tempo
h=L/N;
tau=T/K;
% Inizializzazione del vettore t
t=linspace(0,T,K+1);
% Inizializzazione del vettore x
x=linspace(0,L,N+1);
% Inizializzazione della matrice soluzione u
u=zeros(N+1,K+1);
% Condizione iniziale
u(:,1)=u0(x);
% BC Dirichlet
u(1,:)=c1(t);
% Creation of the matrix A
% We obtained this matrix using the Central Discretization method
e=ones(N-1,1);
A=spdiags([-e,2*e,-e],[-1,0,1],N-1,N-1)/(h^2);
%modify the matrix such that we can use the noimann condition to calculate
%the last node.
%In this case we want to use a second order decentralized approximation
A(end,end-1)=-2/(h^2);
I=speye(N-1,N-1);
for k=1:K
% Assemblaggio termine noto
F=f(x(2:end-1),t(k+1));
% Correzione del termine noto con le condizioni al bordo
F(1)=F(1) + c1(t(k+1))/(h^2);
F(end)=c2*h*2;
% Risoluzione del problema
u(2:end-1,k+1) = ((I+tau*A)\u(2:end-1,k))' + tau*F;
u(end,k+1)=u(end-2,k+1) + c2*2*h;
end

View File

@ -1,107 +0,0 @@
%% Euler forward Dirichlet
L = 1; % domain size (in the lecture notes this is denote h)
T = 2; % time limit (max time)
f = @(x,t) 0*x.*t; % rhs of the more general equation `u_t - d u_xx = f`
c1 = @(t) 1+0*t; % _right_ boundary condition
c2 = @(t) 0*t; % _left_ boundary condition
u0 = @(x) 0*x; % initial values
D = 0.5; % diffusion parameter `d` in `u_t - d u_xx = f`
%uex = @(x,t) cos(x).*exp(t);
N = 10; % nr. of _space_ discretization points
K = 200; % nr. of _time_ discretization points
[x, t, u] = Dirichlet_EA(L, N, T, K, c1, c2, f, u0, D);
% Report stability condition `D Delta T / (Delta x)^2 > 0.5`
Delta_T = T / K;
Delta_x = L / N;
d = D * Delta_T / Delta_x^2;
fprintf("Stability Condition: 0.5 >= D * Delta_T / Delta_x^2 = %f\n", d)
if d > 0.5
fprintf("-> NOT Stable\n")
else
fprintf("-> Stable\n")
end
figure(1)
for ii = 1:K+1 % iterates time
hold on
plot(x, u(:, ii)');
xlim([0 L])
pause(0.05);
hold off
end
% 3D plot of space solution over time
space = linspace(0,L,101);
time = linspace(0,T,201);
[xx,yy] = meshgrid(time,space);
%exsol = uex(yy,xx);
figure(2)
mesh(t,x,u)
%figure(2)
%mesh(xx,yy,exsol)
%% Eulero forward Mixed BC
L=2*pi;
T=5;
f=@(x,t) 0*x.*t;
c1=@(t) 1+0*t;
c2=0;
u0=@(x) 0*x;
D=1.1;
%uex=@(x,t) cos(x).*exp(t);
N=25;
K=200;
[x,t,u]=Mixed_EA(L,N,T,K,c1,c2,f,u0,D);
figure(1)
for ii=1:K+1
plot(x,u(:,ii)');
xlim([0 L])
ylim([0 1.5])
pause(0.02);
end
space=linspace(0,L,101);
time=linspace(0,T,201);
[xx,yy]=meshgrid(time,space);
%exsol=uex(yy,xx);
figure(2)
mesh(t,x,u)
%figure(2)
%mesh(xx,yy,exsol)
%% Eulero Backward Mixed BC
L=2*pi;
T=5;
f=@(x,t) 0*x.*t;
c1=@(t) 1+0*t;
c2=0;
u0=@(x) 0*x;
D=1;
%uex=@(x,t) cos(x).*exp(t);
N=25;
K=200;
[x,t,u]=Mixed_EI(L,N,T,K,c1,c2,f,u0,D);
figure(1)
for ii=1:K+1
plot(x,u(:,ii)');
xlim([0 L])
ylim([0 1.5])
pause(0.02);
end
space=linspace(0,L,101);
time=linspace(0,T,201);
[xx,yy]=meshgrid(time,space);
%exsol=uex(yy,xx);
figure(2)
mesh(t,x,u)
%figure(2)
%mesh(xx,yy,exsol)

View File

@ -1,43 +0,0 @@
%% exercise number 2
%problem characteristics
U=0.7;%convection velocity
delta_x=0.01; %space discretization
L=1; %max length of our domain
T=1; %max time considered
N=L/delta_x; %number of intervals in space
K=100; %number of intervals in time
%two initial conditions
square_pulse=@(x) heaviside(x-0.1) - heaviside(x- 0.3);
gauss_signal=@(x) exp(-10*(4*x - 1).^2);
[x,t,c] = UW_scheme(L,N,T,K,U,square_pulse);
ex_sol=@(x,t) square_pulse(x-U*t);
f1=figure(1);
for ii=1:K+1
clf(f1)
hold on
plot(x,c(:,ii)','-bo');
plot(x,ex_sol(x,t(ii)));
hold off
legend("aproximated solution","exact solution")
xlim([0 L])
ylim([0 1.1])
pause(0.02);
end
[x,t,c] = UW_scheme(L,N,T,K,U,gauss_signal);
ex_sol=@(x,t) gauss_signal(x-U*t);
f2=figure(2);
for ii=1:K+1
clf(f2)
hold on
plot(x,c(:,ii)','-bo');
plot(x,ex_sol(x,t(ii)));
legend("aproximated solution","exact solution")
xlim([0 L])
ylim([0 1.1])
pause(0.02);
end

View File

@ -1,51 +0,0 @@
function [x,t,c] = UW_scheme(xf,N,T,K,U,c0)
% ---- Numerical solution of the linear advection equation in a periodic domain ----
% c_t+ U * c_x = 0 with domain [0,xf]
% given the initial condition c0.
% -----------------------------------------------
% Sintax:
% [x,t,c] = UW_scheme(xf,N,T,K,U,c0)
%
% Input:
% xf end of our domain
% N number of space intervals
% T max time
% K number of time intervals
% U convection velocity
% c0 initial condition
%
% Output:
% x vector of spatial nodes
% t vector of time nodes
% c numerical solution
% Space and time intervals size
dx=xf/N;
dt=T/K;
% initialization of x and t vectors (nodes)
x=linspace(0,xf,N+1)';
t=linspace(0,T,K+1)';
% Solution matrix
c=zeros(N+1,K+1);
% Initial conditions
c(:,1) = c0(x);
% Creating our matrix
e = ones(N+1,1);
B = spdiags([-e,e],[-1,0],N+1,N+1);
I = speye(N+1);
%"printing" courant number and numerical viscosity
C0=(U*dt/dx)
numerical_viscosity= U*dx*(1-C0)/2
%final metrix to compute the solution
A = I - C0*B;
%finding the solution
for k=1:K
c(1:end,k+1) = A*c(1:end,k);
end

View File

@ -1,60 +0,0 @@
function [x,t,u] = calore_CN(L,N,T,K,c1,c2,f,u0)
% ---- Risoluzione dell'equazione del calore ----
% u_t - u_xx = f nell'intervallo [-L,L]
% con condizioni al bordo di Dirichlet
% e condizioni iniziali.
% -----------------------------------------------
% Sintassi:
% [x,t,u]=calore_template(L,N,T,K,c1,c2,fun,u0)
%
% Input:
% L semiampiezza intervallo spaziale (-L,L)
% N numero di sottointervalli in (-L,L)
% T estremo finale intervallo temporale (0,T)
% K numero di sottointervalli in (0,T)
% c1 funzione che descrive la condizione di Dirichlet in x=-L
% c2 funzione che descrive la condizione di Dirichlet in x=L
% f funzione che descrive il termine noto dell'equazione
% u0 funzione che descrive la condizione iniziale in t=0
%
% Output:
% x vettore dei nodi spaziali
% t vettore dei nodi temporali
% u soluzione numerica
% del problema
% Calcolo passo di discretizzazione in spazio e tempo
h=2*L/N;
tau=T/K;
% Inizializzazione del vettore t
t=linspace(0,T,K+1)';
% Inizializzazione del vettore x
x=linspace(-L,L,N+1);
% Inizializzazione della matrice soluzione u
u=zeros(N+1,K+1);
% Condizione iniziale
u(:,1)=u0(x);
% Condizioni al bordo
u(1,:)=c1(t);
u(end,:)=c2(t);
% Costruzione della matrice A
e=ones(N-1,1);
A=spdiags([-e,2*e,-e],[-1,0,1],N-1,N-1)/(h^2);
I=speye(N-1,N-1);
% Ciclo iterativo
for k=1:K
% Assemblaggio termine noto
F1=f(x(2:end-1),t(k));
F2=f(x(2:end-1),t(k+1))
% Correzione del termine noto con le condizioni al bordo
F1(1)=F1(1) + c1(t(k))/(h^2);
F1(end)=F1(end) + c2(t(k))/(h^2);
F2(1)=F2(1) + c1(t(k+1))/(h^2);
F2(end)=F2(end) + c2(t(k+1))/(h^2);
% Risoluzione del problema
u(2:end-1,k+1) = ((I + 0.5*tau*A)\((I - 0.5*tau*A)*u(2:end-1,k)))' + 0.5*tau*F1 + 0.5*tau*F2;
end

BIN
Exercise_03/plots/task01_1_00023.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task01_1_00199.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task01_1_Cinf.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task01_2_00023.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task01_2_00199.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task01_2_Cinf.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task01_3_00023.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task01_3_00199.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task01_4_00023.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task01_4_00199.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task01_all.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task01_x_00023.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task01_x_00199.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task02_gauss_C0_0-7_00150.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task02_gauss_C0_1-1_00150.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task02_gauss_C0_1_00150.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task02_square_C0_0-7_00150.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task02_square_C0_1-1_00150.png (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/plots/task02_square_C0_1_00150.png (Stored with Git LFS) Normal file

Binary file not shown.

View File

@ -12,7 +12,7 @@
% Document meta into % Document meta into
\title{NSSC 2 - Assignement 3} \title{NSSC 2 - Assignement 3}
\author{Bianchi Riccardo, Kapla Daniel, Kuen Jakob, Müller David} \author{Bianchi Riccardo, Kapla Daniel, Kuen Jakob, Müller David}
\date{November 24, 2021} \date{May 17, 2022}
% Set PDF title, author and creator. % Set PDF title, author and creator.
\AtBeginDocument{ \AtBeginDocument{
\hypersetup{ \hypersetup{
@ -88,7 +88,7 @@ We are given the 1D unsteady diffusion equation
\end{equation} \end{equation}
with a scalar diffusion coefficient $D = 10^{-6}$. The domain size is $h$ which is discretized with $N_x$ points to give the grid space $\Delta x$ with a time step $\Delta t$. The initial condition is $C = 0$ inside the domain. with a scalar diffusion coefficient $D = 10^{-6}$. The domain size is $h$ which is discretized with $N_x$ points to give the grid space $\Delta x$ with a time step $\Delta t$. The initial condition is $C = 0$ inside the domain.
Furthermore, denote with $N_t\in\mathbb{N}$ the number of time steps of the simulation where $t_0 = 0$ is the initial time. All $N_t$ discretized time points $t_n = n \Delta t$ leading to the final time point of the simulation as $t_{N_t} = N_t \Delta t$. The discretized $x$ points of the domain are $x_i = (i - 1)\Delta x = (i - 1)h / (N_x - 1)$ for $i = 1, ..., N_x$. This gives exactly $N_x$ grid points with equal distance between adjacent points on the domain $[0, h]$ with the first point $x_1 = 0$ and the last point $x_{N_x} = h$. Furthermore, denote with $N_t\in\mathbb{N}$ are the number of time steps of the simulation where $t_0 = 0$ is the initial time. All $N_t$ discretized time points $t_n = n \Delta t$ leading to the final time point of the simulation as $t_{N_t} = N_t \Delta t$. The discretized $x$ points of the domain are $x_i = (i - 1)\Delta x = (i - 1)h / (N_x - 1)$ for $i = 1, ..., N_x$. This gives exactly $N_x$ grid points with equal distance between adjacent points on the domain $[0, h]$ with the first point $x_1 = 0$ and the last point $x_{N_x} = h$.
In the following we will also use the short hand notation In the following we will also use the short hand notation
\begin{displaymath} \begin{displaymath}
@ -104,7 +104,7 @@ Furthermore, define
\end{displaymath} \end{displaymath}
\subsection{Explicit scheme with Dirichlet/Neumann BC}\label{sec:task01_1} \subsection{Explicit scheme with Dirichlet/Neumann BC}\label{sec:task01_1}
Given the Dirichlet boundary conditions $C(0, t) = 0$ and the Neumann BC $\partial_x C(h, t) = 0$. First we derive a $2^{nd}$ order \emph{explicit} finite difference approach with a $1^{st}$ order discretization in time (see Figure~\ref{fig:ex1}). Given the Dirichlet boundary conditions $C(0, t) = 0$ and the Neumann boundary conditions $\partial_x C(h, t) = 0$. First we derive a $2^{nd}$ order \emph{explicit} finite difference approach with a $1^{st}$ order discretization in time (see Figure~\ref{fig:ex1}).
\begin{figure}[h!] \begin{figure}[h!]
\centering \centering
@ -161,7 +161,7 @@ To derive a second order scheme (assuming $C$ is three time continuously differe
\begin{displaymath} \begin{displaymath}
C(x + \Delta x, t) = C(x, t) + \partial_x C(x, t)\Delta x + \partial_x^2 C(x, t)\frac{\Delta x^2}{2} + \partial_x^3 C(x, t)\frac{\Delta x^3}{6} + \mathcal{O}(\Delta x^4). C(x + \Delta x, t) = C(x, t) + \partial_x C(x, t)\Delta x + \partial_x^2 C(x, t)\frac{\Delta x^2}{2} + \partial_x^3 C(x, t)\frac{\Delta x^3}{6} + \mathcal{O}(\Delta x^4).
\end{displaymath} \end{displaymath}
Replace $\Delta x$ with $\pm\Delta x$ and add the two equation together, then Replace $\Delta x$ with $\pm\Delta x$ and add the two equation together which gives
\begin{displaymath} \begin{displaymath}
C(x + \Delta x, t) + C(x - \Delta x, t) = 2 C(x, t) + \partial_x^2 C(x, t)\Delta x^2 + \mathcal{O}(\Delta x^4). C(x + \Delta x, t) + C(x - \Delta x, t) = 2 C(x, t) + \partial_x^2 C(x, t)\Delta x^2 + \mathcal{O}(\Delta x^4).
\end{displaymath} \end{displaymath}
@ -171,7 +171,7 @@ The first and third order terms drop out due to different signs. Finally, this r
\end{displaymath} \end{displaymath}
which is an approximation of the second derivative with error proportional to $\Delta x^2$. For the first order approximation in time the same trick can be used except that only the first two terms of the Taylor expantion need to be considured. which is an approximation of the second derivative with error proportional to $\Delta x^2$. For the first order approximation in time the same trick can be used except that only the first two terms of the Taylor expantion need to be considured.
Therefore, the $2^{nd}$ order scheme for space and $1^{st}$ order in time at $(i, n)$ using the explicit scheme derives as Therefore, the $2^{nd}$ order scheme for space and $1^{st}$ order in time at $(i, n)$ using the explicit scheme derive as
\begin{displaymath} \begin{displaymath}
\partial_x^2 \widehat{C}_i^n = \frac{\widehat{C}_{i-1}^{n} - 2\widehat{C}_{i}^{n} + \widehat{C}_{i+1}^{n}}{\Delta x^2}, \partial_x^2 \widehat{C}_i^n = \frac{\widehat{C}_{i-1}^{n} - 2\widehat{C}_{i}^{n} + \widehat{C}_{i+1}^{n}}{\Delta x^2},
\qquad \qquad
@ -181,7 +181,7 @@ Substitution into \eqref{eq:ex1} yields after rearranging the update rule for in
\begin{equation}\label{eq:task01_1_update} \begin{equation}\label{eq:task01_1_update}
\widehat{C}_i^{n+1} = \widehat{C}_i^n + \frac{D\Delta t}{\Delta x^2}(\widehat{C}_{i-1}^{n} - 2\widehat{C}_{i}^{n} + \widehat{C}_{i+1}^{n}). \widehat{C}_i^{n+1} = \widehat{C}_i^n + \frac{D\Delta t}{\Delta x^2}(\widehat{C}_{i-1}^{n} - 2\widehat{C}_{i}^{n} + \widehat{C}_{i+1}^{n}).
\end{equation} \end{equation}
This holds at $x_i$ where $i = 2, ..., N_x - 1$, or in other words everywhere inside the space domain excluding the boundary. The boundary needs to be handled seperately. The left boundary condition is a Dirichlet constraint which is simply constanct giving $C_1^n = 1$ for all time while the Neumann condition on the right requires an additional discretization step for computing the next value. To achieve a $2^{nd}$ order discretization scheme it's usualy required to take the sum/difference of an Taylor expantion into the left and right direction at a given point. As we are on the boundary there is no point at the right, which can be simply solved by adding a ghost point $x_{N_x+1}$ to the system (which will not show up in the final scheme). This leads to the two $3^{rd}$ order expantions evaluated at the right boundary as This holds at $x_i$ where $i = 2, ..., N_x - 1$, or in other words, everywhere inside the space domain except the boundary. The boundary needs to be handled seperately. The left boundary condition is a Dirichlet constraint which is simply a constant giving $C_1^n = 1$ for all time while the Neumann condition on the right requires an additional discretization step for computing the next value. To achieve a $2^{nd}$ order discretization scheme it is usualy required to take the sum/difference of an Taylor expantion to the left and right direction at a given point.However, on the boundary there is no point at the right, which can be simply solved by adding a ghost point $x_{N_x+1}$ to the system (which will not show up in the final scheme). This leads to the two $3^{rd}$ order expantions evaluated at the right boundary as
\begin{align*} \begin{align*}
C_{N_x-1}^n &= C_{N_x}^n - \partial_x C_{N_x}^n\Delta x + \partial_x^2 C_{N_x}^n + \mathcal{O}(\Delta x^3), \\ C_{N_x-1}^n &= C_{N_x}^n - \partial_x C_{N_x}^n\Delta x + \partial_x^2 C_{N_x}^n + \mathcal{O}(\Delta x^3), \\
C_{N_x+1}^n &= C_{N_x}^n + \partial_x C_{N_x}^n\Delta x + \partial_x^2 C_{N_x}^n + \mathcal{O}(\Delta x^3). C_{N_x+1}^n &= C_{N_x}^n + \partial_x C_{N_x}^n\Delta x + \partial_x^2 C_{N_x}^n + \mathcal{O}(\Delta x^3).
@ -194,7 +194,7 @@ which yields a $2^{nd}$ order equation for the Neumann BC on the right as
\begin{displaymath} \begin{displaymath}
C_{N_x-1}^n = C_{N_x+1}^n. C_{N_x-1}^n = C_{N_x+1}^n.
\end{displaymath} \end{displaymath}
Now where ready to give an equation system of the entire system (including the boundaryies) which reads This results in an system of equations containing the entire system (including the boundaryies) which reads
\begin{align*} \begin{align*}
\widehat{C}_1^{n+1} &= 1, \\ \widehat{C}_1^{n+1} &= 1, \\
\widehat{C}_i^{n+1} &= d\widehat{C}_{i-1}^{n} + (1 - 2d)\widehat{C}_{i}^{n} + d\widehat{C}_{i+1}^{n}, \widehat{C}_i^{n+1} &= d\widehat{C}_{i-1}^{n} + (1 - 2d)\widehat{C}_{i}^{n} + d\widehat{C}_{i+1}^{n},
@ -218,10 +218,19 @@ for a tridiagonal matrix $A_1$ of dimensions $N_x\times N_x$
\end{pmatrix}. \end{pmatrix}.
\end{displaymath} \end{displaymath}
\subsubsection{Stability} \todo{Prove that the update scheme is unstable for $d = \frac{D\Delta x^2}{\Delta t} > 0.5$} \begin{figure}[!h]
\centering
\includegraphics[width=0.48\textwidth]{task01_1_00023.png}
\includegraphics[width=0.48\textwidth]{task01_x_00023.png}
\includegraphics[width=0.48\textwidth]{task01_1_00199.png}
\includegraphics[width=0.48\textwidth]{task01_x_00199.png}
\caption{\label{fig:task01_stable}Comparison for stable (left, $d = 0.48 \leq 0.5$) vs. unstable (right, $d = 0.55 > 0.5$). At the top the simulation snapshot after 120 time steps and at the bottom after 1000 timesteps.}
\end{figure}
\subsection{Explicit scheme with Dirichlet BC at both boundaries}\label{sec:task01_2} \subsection{Explicit scheme with Dirichlet BC at both boundaries}\label{sec:task01_2}
This is identical to the previous task except the right boundary condition. The new boundary condition is a Dirichlet BC with $C_{N_x}^n = 0$ which alters the equation system from Section~\ref{sec:task01_1} in the last two equations to In this task two Dirichlet boundary conditions are given. The new boundary condition is a Dirichlet BC with $C_{N_x}^n = 0$ which alters the equation system from Section~\ref{sec:task01_1} in the last two equations to
\begin{align*} \begin{align*}
\widehat{C}_1^{n+1} &= 1, \\ \widehat{C}_1^{n+1} &= 1, \\
\widehat{C}_i^{n+1} &= d\widehat{C}_{i-1}^{n} + (1 - 2d)\widehat{C}_{i}^{n} + d\widehat{C}_{i+1}^{n}, \widehat{C}_i^{n+1} &= d\widehat{C}_{i-1}^{n} + (1 - 2d)\widehat{C}_{i}^{n} + d\widehat{C}_{i+1}^{n},
@ -229,7 +238,7 @@ This is identical to the previous task except the right boundary condition. The
\widehat{C}_{N_x-1}^{n+1} &= d\widehat{C}_{N_x-2}^{n} + (1 - 2d)\widehat{C}_{N_x-1}^{n}, \\ \widehat{C}_{N_x-1}^{n+1} &= d\widehat{C}_{N_x-2}^{n} + (1 - 2d)\widehat{C}_{N_x-1}^{n}, \\
\widehat{C}_{N_x}^{n+1} &= 0. \widehat{C}_{N_x}^{n+1} &= 0.
\end{align*} \end{align*}
which yields a new tridiagonal matrix $A_2$ of size $N_x\times N_x$ for the update rule $C^{n+1} = A_2 C^{n}$ as which yields to a new tridiagonal matrix $A_2$ of size $N_x\times N_x$ for the update rule $C^{n+1} = A_2 C^{n}$ as
\begin{displaymath} \begin{displaymath}
A_2 = \begin{pmatrix} A_2 = \begin{pmatrix}
1 & 0 \\ 1 & 0 \\
@ -242,10 +251,17 @@ which yields a new tridiagonal matrix $A_2$ of size $N_x\times N_x$ for the upda
\end{pmatrix}. \end{pmatrix}.
\end{displaymath} \end{displaymath}
\todo{include pictures and compare with Section~\ref{sec:task01_1}} See Figure~\ref{fig:task01_inf} for $\lim_{t\to\infty}C(x, t)$ simulation (many timesteps) for subtask 1 and 2. In comparison both lead to a straight line solution (last energy) such that the Dirichlet or Neumann BC are fullfilled.
\begin{figure}[!h]
\centering
\includegraphics[width=0.48\textwidth]{task01_1_Cinf.png}
\includegraphics[width=0.48\textwidth]{task01_2_Cinf.png}
\caption{\label{fig:task01_inf}Solution from Section~\ref{sec:task01_1} (left) and Section~\ref{sec:task01_2} (right) for $t\to\infty$ simulated by setting number of time steps to $N_t = 20000$.}
\end{figure}
\subsection{Implicit scheme with Dirichlet/Neumann BC}\label{sec:task01_3} \subsection{Implicit scheme with Dirichlet/Neumann BC}\label{sec:task01_3}
The following computations are similar in nature to the explicit scheme, therefore well keep it short. The implicit discretization for the derivatives is (see: Figure~\ref{fig:ex1}) then The following computations are similar in nature to the explicit scheme. The implicit discretization for the derivatives are (see: Figure~\ref{fig:ex1})
\begin{displaymath} \begin{displaymath}
\partial_x^2 \widehat{C}_{i}^{n+1} = \frac{\widehat{C}_{i-1}^{n+1} - 2\widehat{C}_{i}^{n+1} + \widehat{C}_{i+1}^{n+1}}{\Delta x^2}, \partial_x^2 \widehat{C}_{i}^{n+1} = \frac{\widehat{C}_{i-1}^{n+1} - 2\widehat{C}_{i}^{n+1} + \widehat{C}_{i+1}^{n+1}}{\Delta x^2},
\qquad \qquad
@ -275,9 +291,13 @@ where the $N_x\times N_x$ tridiagonal matrix $A_3$ has the form
\end{displaymath} \end{displaymath}
and the update is performed by solving for $C^{n+1}$. and the update is performed by solving for $C^{n+1}$.
\todo{include pictures and compare with Section~\ref{sec:task01_1}} \begin{figure}[!h]
\centering
\includegraphics[width = 0.8\textwidth]{task01_all.png}
\caption{\label{fig:task01_all}Comparison of all four subtasks. As expected 1.1, 1.3 and 1.4 solve the same probelm and as such have (basically) identical solution while the subtask 1.2 solves for different BC (Dirichlet) which is fullfilled.}
\end{figure}
\subsection{Second order in time for implicit scheme with Dirichlet/Neumann BC} \subsection{Second order in time for implicit scheme with Dirichlet/Neumann BC}\label{sec:task01_4}
\begin{figure}[h!] \begin{figure}[h!]
\centering \centering
\begin{tikzpicture}[>=latex] \begin{tikzpicture}[>=latex]
@ -302,12 +322,12 @@ and the update is performed by solving for $C^{n+1}$.
minimum size = 4pt] at (\x, \y) {}; minimum size = 4pt] at (\x, \y) {};
} }
\end{tikzpicture} \end{tikzpicture}
\caption{\label{fig:ex4}.Crank-Nicolson dependency relations.} \caption{\label{fig:ex4}Crank-Nicolson dependency relations.}
\end{figure} \end{figure}
To derive the second order accurate scheme in time we employ the Crank-Nicolson approach to derive the descritization at the time midpoints $t_{n+1/2} = (n + 1/2)\Delta t$ which we index with $n+1/2$. To derive the second order accurate scheme in time one uses the Crank-Nicolson approach to derive the descritization at the time midpoints $t_{n+1/2} = (n + 1/2)\Delta t$ which are indexed with $n+1/2$.
By expantion into $\pm\Delta t / 2$ at the midpoints we get By expantion into $\pm\Delta t / 2$ at the midpoints one get
\begin{align*} \begin{align*}
C_i^n &= C_i^{n+1/2} - \partial_tC_i^{n+1/2}\frac{\Delta t}{2} + \partial_t^2 C_i^{n+1/2}\frac{\Delta t^2}{8} + \mathcal{O}(\Delta t^3), \\ C_i^n &= C_i^{n+1/2} - \partial_tC_i^{n+1/2}\frac{\Delta t}{2} + \partial_t^2 C_i^{n+1/2}\frac{\Delta t^2}{8} + \mathcal{O}(\Delta t^3), \\
C_i^{n+1} &= C_i^{n+1/2} + \partial_tC_i^{n+1/2}\frac{\Delta t}{2} + \partial_t^2 C_i^{n+1/2}\frac{\Delta t^2}{8} + \mathcal{O}(\Delta t^3). C_i^{n+1} &= C_i^{n+1/2} + \partial_tC_i^{n+1/2}\frac{\Delta t}{2} + \partial_t^2 C_i^{n+1/2}\frac{\Delta t^2}{8} + \mathcal{O}(\Delta t^3).
@ -325,12 +345,12 @@ leading to an second order in time and space discretization scheme
\partial_t \widehat{C}_{i}^{n+1/2} - D \partial_x^2 \widehat{C}_{i}^{n+1/2} = \mathcal{O}(\Delta x^2 + \Delta t^2). \partial_t \widehat{C}_{i}^{n+1/2} - D \partial_x^2 \widehat{C}_{i}^{n+1/2} = \mathcal{O}(\Delta x^2 + \Delta t^2).
\end{displaymath} \end{displaymath}
Next we need to ensure the evaluation to be done at the known grid points instead of the unknown time midpoints. Therefore, we replace the evaluation at the time midpoints by there in time mean of the two direct neighbouring grid points Next one need to ensure the evaluation to be done at the known grid points instead of the unknown time midpoints. Therefore, replacing the evaluation at the time midpoints by the time mean of the two direct neighbouring grid points
\begin{displaymath} \begin{displaymath}
\widehat{C}_i^{n+1/2} = \frac{\widehat{C}_i^{n} + \widehat{C}_i^{n+1}}{2}. \widehat{C}_i^{n+1/2} = \frac{\widehat{C}_i^{n} + \widehat{C}_i^{n+1}}{2}.
\end{displaymath} \end{displaymath}
After a bit of algebra and condideration of the Dirichlet and Neumann BC we end up with After a bit of algebra and consideration of the Dirichlet and Neumann BC one end up with
\begin{align*} \begin{align*}
\widehat{C}_1^{n+1} &= \widehat{C}_1^{n} = 1, \\ \widehat{C}_1^{n+1} &= \widehat{C}_1^{n} = 1, \\
\widehat{C}_i^{n+1} - \frac{d}{2}\widehat{C}_{i-1}^{n+1} + dC_{i}^{n+1} - \frac{d}{2}\widehat{C}_{i+1}^{n+1} \widehat{C}_i^{n+1} - \frac{d}{2}\widehat{C}_{i-1}^{n+1} + dC_{i}^{n+1} - \frac{d}{2}\widehat{C}_{i+1}^{n+1}
@ -339,7 +359,7 @@ After a bit of algebra and condideration of the Dirichlet and Neumann BC we end
\widehat{C}_i^{n+1} - d\widehat{C}_{N_x-1}^{n+1} + dC_{N_x}^{n+1} \widehat{C}_i^{n+1} - d\widehat{C}_{N_x-1}^{n+1} + dC_{N_x}^{n+1}
&= C_i^{n} + d\widehat{C}_{N_x-1}^{n} - dC_{N_x}^{n}. &= C_i^{n} + d\widehat{C}_{N_x-1}^{n} - dC_{N_x}^{n}.
\end{align*} \end{align*}
To write this again in a compact matrix equation system we define the $N_x\times N_x$ matrix $A_4$ as To write this again in a compact matrix equation system one defines the $N_x\times N_x$ matrix $A_4$ as
\begin{displaymath} \begin{displaymath}
A_4 = \begin{pmatrix} A_4 = \begin{pmatrix}
0 & 0 \\ 0 & 0 \\
@ -356,4 +376,79 @@ then the entire system has the form
\end{displaymath} \end{displaymath}
where $I$ is the identity. Again, the update is performed by solving this tridiagonal system for $C^{n+1}$. where $I$ is the identity. Again, the update is performed by solving this tridiagonal system for $C^{n+1}$.
For comparison against the other methods see Figure~\ref{fig:task01_all}.
\subsection{Notes on Implementation}
All four tasks follow the same basic sructure which can be written as a generalized linear system of the form
\begin{displaymath}
L C^{n+1} = R C^n
\end{displaymath}
for updating the solution by one time step. The two matrices are dependent on the used scheme and the boundary conditions, see Table~\ref{tab:ex1}.
\begin{table}[h!]
\centering
\begin{tabular}{l | c c}
Task & $L$ & $R$ \\ \hline
1.1 & $I$ & $A_1$ \\
1.2 & $I$ & $A_2$ \\
1.3 & $A_3$ & $I$ \\
1.4 & $I + A_4$ & $I - A_4$ \\
\end{tabular}
\caption{\label{tab:ex1}Left and right hand side matrices of the general update rule $L C^{n+1} = R C^n$ with $A_i$ depending on the subtask as defined in Section~\ref{sec:task01_1} till \ref{sec:task01_4} and $I$ the identity matrix.}
\end{table}
As a boild down solver the following is sufficient to compute the solution given an initial solution $C_0$ (argument \texttt{C}) with the appropriately defined matrices $L, R$ (arguments \texttt{L, R}) which incorporate the scheme, space and time step size $\Delta x, \Delta t$ as well as the dispersion parameter $D$ and the boundary conditions) at time $T = N_t \Delta$ (implicitly defined through the parameter \texttt{nx}).
\begin{python}
def solve(C, L, R, nt):
for t in range(nt):
rhs = C if R is None else R @ C
C = rhs if L is None else solve_banded((1, 1), L.data, rhs)
return C
\end{python}
This code assumes \texttt{C} to be a 1D \texttt{numpy} array as well as $L, R$ to by \texttt{scipy.sparse.spmatrix} sparse tridiagonal matrices or \texttt{None} which is equivalent to the identity matrix.
The additional code in the script \texttt{task01.py} contains setting up the needed matrices and plotting.
\newpage
\section{Numerical solution of the linear advection equation in a periodic domain}
First a few words about the upwind scheme. With a positive $U$ the upwind means that it is discretized by ``going left'', meaning into the opposite direction of the current position. In a similar faschion as above this leads to the discretization
\begin{align*}
\partial_x C_i^{n} &= \frac{C_i^n - C_{i-1}^n}{\Delta x} + \mathcal{O}(\Delta x), \\
\partial_t C_i^{n+1} &= \frac{C_i^{n+1} - C_i^n}{\Delta t} + \mathcal{O}(\Delta t).
\end{align*}
One ends up with a first order discretized equation
\begin{displaymath}
0 = \frac{\widehat{C}_i^{n+1} - \widehat{C}_i^n}{\Delta t} + U \frac{\widehat{C}_i^n - \widehat{C}_{i-1}^n}{\Delta x}
\end{displaymath}
leading to the update rule
\begin{align*}
C_i^{n+1} &= (1 - C_0)C_i^n + C_oC_{i-1}^n, && i = 2, ..., N_x
\end{align*}
with the current $C_0 = \frac{U \Delta x}{\Delta t}$. The boundary condition enforces periodizity which means $C_1^n = C_{N_x}^n$ and therefore the update for the left most point is
\begin{displaymath}
C_1^{n+1} = (1 - C_0)C_1^n + C_oC_{N_x}^n.
\end{displaymath}
\begin{figure}[!h]
\centering
\includegraphics[width = 0.3\textwidth]{task02_gauss_C0_0-7_00150.png}
\includegraphics[width = 0.3\textwidth]{task02_gauss_C0_1_00150.png}
\includegraphics[width = 0.3\textwidth]{task02_gauss_C0_1-1_00150.png}
\includegraphics[width = 0.3\textwidth]{task02_square_C0_0-7_00150.png}
\includegraphics[width = 0.3\textwidth]{task02_square_C0_1_00150.png}
\includegraphics[width = 0.3\textwidth]{task02_square_C0_1-1_00150.png}
\caption{\label{fig:task02_all}For initial Gauss wavelet (top) and Square wavelet (bottom) after 150 iterations (one leap around the boundary) with different current $C_0$. Left: $C_0 = 0.7$, Center: $C_0 = 1$ (perfect) and Right: $C_0 = 1.1$ (unstable).}
\end{figure}
For implementation one gets a very simple solver which computes the wave form at time $T = N_t \Delta t$ with $\Delta t = U \Delta x / C_0$ as the following;
\begin{python}
def solve(C, C0, nt):
nx = C.shape[0]
im1 = np.array([nx - 1] + list(range(nx - 1)))
for t in range(nt):
C = (1 - C0) * C + C0 * C[im1]
return C
\end{python}
where \texttt{C} is the initial wave form, the constant \texttt{C0} is the current $C_0$ and \texttt{nt} is the number of time steps to be performed. Again, the attached script \texttt{task02.py} implements this with additional code for setting up the initial conditions and plotting.
\end{document} \end{document}

154
Exercise_03/task01.py Normal file
View File

@ -0,0 +1,154 @@
import numpy as np
from scipy.sparse import spdiags
from scipy.linalg import solve_banded
from matplotlib import pyplot as plt
# Config
D = 1e-6 # diffusion coefficient
h = 1 # space domain (max x size)
T = 2e5 # solution end time
nx = 50 # nr of space discretization points
nt = 1000 # nr of time discretization points
# derived constants
dx = h / (nx - 1) # space step size
dt = T / nt # time step size
d = dt * D / dx**2 # stability/stepsize coefficient
# Task 1.1 equation system matrix `C^n+1 = A1 C^n`
def A1(d):
a1 = np.zeros([3, nx])
# superdiagonal (above the main diagonal)
a1[0, 1] = 0
a1[0, 2:] = d
# main diagonal
a1[1, 0] = 1
a1[1, 1:] = 1 - 2 * d
# subdiagonal (below the main diagonal)
a1[2, :-2] = d
a1[2, -2] = 2 * d
# convert to sparse tridiagonal matrix
return spdiags(a1, [1, 0, -1], nx, nx)
# Task 1.2 equation system matrix `C^n+1 = A2 C^n`
def A2(d):
a2 = A1(d).data
# addapt boundary condition coefficients, the rest is identical
a2[0, -1] = a2[1, -1] = a2[2, -2] = 0
# convert to sparse tridiag
return spdiags(a2, [1, 0, -1], nx, nx)
# Task 1.3 equation system matrix for update rule `A3 C^n+1 = C^n`
def A3(d):
a3 = np.zeros([3, nx])
# superdiagonal
a3[0, 1] = 0
a3[0, 2:] = -d
# main diagonal
a3[1, 0] = 1
a3[1, 1:] = 1 + 2 * d
# supdiagonal
a3[2, :-2] = -d
a3[2, -2] = -2 * d
# as tridiag
return spdiags(a3, [1, 0, -1], nx, nx)
# Task 1.4 equation system matrix `A4`` in `(I + A4) C^n+1 = (I - A4) C^n`
def A4(d):
a4 = np.zeros([3, nx])
# superdiagonal
a4[0, 2:] = -d / 2
# main diagonal
a4[1, 1:] = d
# supdiagonal
a4[2, :-2] = -d / 2
a4[2, -2] = -d
return spdiags(a4, [1, 0, -1], nx, nx)
# setup `I + A4` and `I - A4`
def IpA4(d):
ipA4 = A4(d).data
ipA4[1, ] += 1
return spdiags(ipA4, [1, 0, -1], nx, nx)
def ImA4(d):
imA4 = -A4(d).data
imA4[1, ] += 1
return spdiags(imA4, [1, 0, -1], nx, nx)
# plots a solution to file, simply for code size reduction in `solve`
def plot(x, C, name):
plt.clf()
plt.plot(x, C)
plt.xlim([0, h])
plt.ylim([0, 1.2])
plt.title(name)
plt.xlabel("Distance from Source: x")
plt.ylabel("Concentration: C(x, t)")
plt.savefig(f"plots/{name}_{plot.i:0>5}.png")
plot.i += 1
plot.i = 0
# solver for `C_t - D C_xx = 0` given discretized equation system matrices for
# a time update rule `L C^n+1 = R C^n` starting from an initial solution `C`.
# The update is iterated for `nt` time steps corresponding to a solution at
# time `T = dt nt`
def solve(C, L, R, nt, plotname = None, plotskip = 5):
if plotname is not None:
plot.i = 0 # (re)set plot counter
x = np.linspace(0, h, nx) # space grid points
fig = plt.figure(figsize = (8, 6), dpi = 100) # setup plot figure size
plot(x, C, plotname) # plot initial solution
for t in range(nt):
rhs = C if R is None else R @ C
C = rhs if L is None else solve_banded((1, 1), L.data, rhs)
if plotname is not None and ((t + 1) % plotskip == 0):
plot(x, C, plotname) # plot current solution
return C
# setup initial solution
C0 = np.zeros(nx)
C0[0] = 1
# run simulation and generate plots for each subtask
C1 = solve(C0, None, A1(d), nt, "task01_1")
Cx = solve(C0, None, A1(0.55), 200, "task01_x", 1) # run with unstable `d > 0.5`
C2 = solve(C0, None, A2(d), nt, "task01_2")
C3 = solve(C0, A3(d), None, nt, "task01_3")
C4 = solve(C0, IpA4(d), ImA4(d), nt, "task01_4")
# to convert generated image sequence to video use:
# $> ffmpeg -y -r 60 -i plots/task01_1_%05d.png -pix_fmt yuv420p video_1_1.mp4
# $> ffmpeg -y -r 60 -i plots/task01_x_%05d.png -pix_fmt yuv420p video_1_x.mp4
# $> ffmpeg -y -r 60 -i plots/task01_2_%05d.png -pix_fmt yuv420p video_1_2.mp4
# $> ffmpeg -y -r 60 -i plots/task01_3_%05d.png -pix_fmt yuv420p video_1_3.mp4
# $> ffmpeg -y -r 60 -i plots/task01_4_%05d.png -pix_fmt yuv420p video_1_4.mp4
# compute till no time change, e.g. `t -> inf`
C1inf = solve(C0, None, A1(d), 20000)
C2inf = solve(C0, None, A2(d), 20000)
x = np.linspace(0, h, nx) # space grid points
fig = plt.figure(figsize = (8, 6), dpi = 100) # setup plot figure size
plot.i = 0 # reset plot counter
plot(x, C1inf, "Cinf")
plot(x, C2inf, "Cinf")
# # All plots in one
C1 = solve(C0, None, A1(d), nt)
C2 = solve(C0, None, A2(d), nt)
C3 = solve(C0, A3(d), None, nt)
C4 = solve(C0, IpA4(d), ImA4(d), nt)
x = np.linspace(0, h, nx) # space grid points
plt.figure(figsize = (8, 6), dpi = 100) # setup plot figure size
plt.plot(x, C1, label = "1.1")
plt.plot(x, C2, label = "1.2")
plt.plot(x, C3, label = "1.3")
plt.plot(x, C4, label = "1.4")
plt.xlim([0, h])
plt.ylim([0, 1.2])
plt.title("All")
plt.legend()
plt.xlabel("Distance from Source: x")
plt.ylabel("Concentration: C(x, t)")
plt.savefig(f"plots/task01_all.png")

View File

@ -1,63 +0,0 @@
# Task 1.1, 1.2
import numpy as np
from typing import Callable
from matplotlib import pyplot as plt
# Config
D = 1e-6 # diffusion coefficient
h = 1 # space domain (max x size)
T = 2e6 # solution end time
nx = 50 # nr of space discretization points
nt = 20000 # nr of time discretization points
# derived constants
dx = h / (nx - 1) # space step size
dt = T / (nt - 1) # time step size
d = dt * D / dx**2 # stability/stepsize coefficient
# report stability
if d > 0.5:
print("NOT Stable")
else:
print("Stable")
# explicit scheme integration for `u_t - D u_xx = 0` with boundary conditions
# enforced by `set_bounds` and initial conditions `initial`.
def integrate(*, name: str, initial: np.array, set_bounds: Callable[[np.array], None]) -> None:
C = initial
# Setup boundary conditions
set_bounds(C)
i = 0 # index for plot generation
plt.figure(figsize = (8, 6), dpi = 100)
for t in range(nt):
# every 400'th time step save a plot
if t % (nt // 400) == 0:
plt.clf()
plt.plot(np.linspace(0, h, nx), C)
plt.xlim([0, h])
plt.ylim([0, 1.2])
plt.savefig(f"plots/{name}_{i:0>5}.png")
i += 1
# update solution using the explicit schema
C[1:-1] += d * (C[2:] - 2 * C[1:-1] + C[:-2])
# update right Neumann BC
set_bounds(C)
# Subtask 1 boundary conditions (Dirichlet and Neumann)
def bounds_1(C):
C[0] = 1
C[-1] = C[-2]
# Subtask 2 boundary conditions (two Dirichlet)
def bounds_2(C):
C[0] = 1
C[-1] = 0
# run simulations
integrate(name = 'task01_1', initial = np.zeros(nx), set_bounds = bounds_1)
integrate(name = 'task01_2', initial = np.zeros(nx), set_bounds = bounds_2)
# to convert generated image sequence to video use:
# $> ffmpeg -r 60 -i plots/task01_1_%05d.png -pix_fmt yuv420p video_1_1.mp4
# $> ffmpeg -r 60 -i plots/task01_2_%05d.png -pix_fmt yuv420p video_1_2.mp4

View File

@ -1,50 +0,0 @@
# Task 1.3
import numpy as np
from typing import Callable
from matplotlib import pyplot as plt
# Config
D = 1e-6 # diffusion coefficient
h = 1 # space domain (max x size)
T = 2e6 # solution end time
nx = 50 # nr of space discretization points
nt = 20000 # nr of time discretization points
# derived constants
dx = h / (nx - 1) # space step size
dt = T / (nt - 1) # time step size
d = dt * D / dx**2 # stability/stepsize coefficient
# Setup implicit scheme equation matrix
T = (1 + 2 * d) * np.eye(nx) - d * np.eye(nx, k = 1) - d * np.eye(nx, k = -1)
# fix boundary condition equations
T[0, 0] = 1 # Left Dirichlet BC
T[0, 1] = 0
T[-1, -2] = 1 # Right Neumann BC
T[-1, -1] = 0
# Set initial solution
C = np.zeros(nx)
C[0] = 1
C[-1] = C[-2] # (0 = 0)
i = 0 # index for plot generation
plt.figure(figsize = (8, 6), dpi = 100)
for t in range(nt):
# every 400'th time step save a plot
if t % (nt // 400) == 0:
plt.clf()
plt.plot(np.linspace(0, h, nx), C)
plt.xlim([0, h])
plt.ylim([0, 1.2])
plt.savefig(f"plots/task01_3_{i:0>5}.png")
i += 1
# update solution using the implicit schema
C = np.linalg.solve(T, C)
# fix BC conditions (theoretically, they are set by the update but for
# stability reasons (numerical) we enforce the correct values)
C[0] = 1
C[-1] = C[-2]
# to convert generated image sequence to video use:
# $> ffmpeg -r 60 -i plots/task01_3_%05d.png -pix_fmt yuv420p video_1_3.mp4

View File

@ -1,74 +0,0 @@
# Task 1.4
import numpy as np
import scipy, scipy.sparse, scipy.linalg
from matplotlib import pyplot as plt
# plotting helper function
def plot(C, x, t, name):
plt.clf()
plt.plot(x, C)
plt.xlim([0, x[-1]])
plt.ylim([0, 1.2])
plt.title(f"{name} - time: {t: >12.0f}")
plt.savefig(f"plots/{name}_{plot.fignr:0>5}.png")
plot.fignr += 1
# set static variable as function attribute
plot.fignr = 1
# and initialize a figure
plt.figure(figsize = (8, 6), dpi = 100)
# Config
D = 1e-6 # diffusion coefficient
h = 1 # space domain (max x size)
T = 1e5 # solution end time
nx = 50 # nr of space discretization points
nt = 2000 # nr of time points
# derived constants
dx = h / (nx - 1) # space step size
dt = T / (nt - 1) # time step size
d = dt * D / dx**2 # stability/stepsize coefficient
# setup matrices for implicit update scheme `A C^n+1 = (2 I - A) C^n = B C^n`
# The matrix `A` is represented as a `3 x (nx - 1)` matrix as expected by
# `scipy.linalg.solve_banded` representing the 3 diagonals of the tridiagonal
# matrix `A`.
A = np.zeros((3, nx - 1))
# upper minor diagonal
A[0, 2:] = -d
# main diagonal
A[1, 0] = 1
A[1, 1:] = 1 + 2 * d
# and lower minor diagonal
A[2, :-2] = -d
A[2, -2] = -2 * d
# The right hand side matrix `B = 2 I - A` which is represented as a sparse
# matrix in diag. layout which allows for fast matrix vector multiplication
B = -A
B[1, ] += 2
B = scipy.sparse.spdiags(B, (1, 0, -1), B.shape[1], B.shape[1])
# set descritized space evaluation points `x` of `C`
x = np.linspace(0, h, nx)
# Set initial solution
C = np.zeros(nx)
C[0] = 1
C[-1] = C[-3] # (0 = 0)
for n in range(nt):
# generate roughly 400 plots
if n % (nt // 400) == 0:
plot(C, x, n * dt, "task01_4")
# update solution using the implicit schema
C[:-1] = scipy.linalg.solve_banded((1, 1), A, B @ C[:-1])
# set/enforce boundary conditions
C[0] = 1 # left Dirichlet (theoretically not needed)
C[-1] = C[-3] # right Neumann condition
# plot final state
plot(C, x, T, "task01_4")
# to convert generated image sequence to video use:
# $> ffmpeg -r 60 -i plots/task01_4_%05d.png -pix_fmt yuv420p video_1_4.mp4

68
Exercise_03/task02.py Normal file
View File

@ -0,0 +1,68 @@
import numpy as np
from scipy.sparse import spdiags
from scipy.linalg import solve_banded
from matplotlib import pyplot as plt
# Config
dx = 0.01 # space step size
h = 1.5 # space domain (max x size)
# plots a solution to file, simply for code size reduction in `solve`
def plot(x, t, C, C0, analytic, plotname):
plt.clf()
plt.plot(x, C, "o")
plt.xlim([0, h])
plt.title(plotname)
plt.xlabel("Periodic Space: x")
plt.ylabel("Wave Form: C(x, t)")
if analytic is not None:
plt.plot(x, analytic(x, t, C0))
plt.savefig(f"plots/{plotname}_{plot.i:0>5}.png")
plot.i += 1
plot.i = 0
# solver for `C_t + U C_x = 0` given initial state `C`, the current `C0` which
# incorporates `U` as well as the space and time discretization and the number
# of time steps `nt`.
def solve(C, C0, nt, plotname = None, plotskip = 1, analytic = None):
nx = C.shape[0]
if plotname is not None:
plot.i = 0 # (re)set plot counter
x = np.linspace(0, h, nx) # space grid points
fig = plt.figure(figsize = (8, 6), dpi = 100) # setup plot figure size
plot(x, 0, C, C0, analytic, plotname) # plot initial solution
# space grid point `Indices Minus 1`
im1 = np.array([nx - 1] + list(range(nx - 1)))
for t in range(nt):
# upwind update scheme
C = (1 - C0) * C + C0 * C[im1]
if plotname is not None and ((t + 1) % plotskip == 0):
plot(x, t + 1, C, C0, analytic, plotname) # plot current solution
return C
# setup initial conditions
x = np.linspace(0, h, int(h / dx))
gauss = np.exp(-10 * (4 * x - 1)**2)
square = np.array((0.1 < x) * (x < 0.3), dtype = "float")
def gauss_analytic(x, t, C0):
return np.exp(-10 * (4 * (np.mod(x - t * dx / C0, h + 1e-9)) - 1)**2)
def square_analytic(x, t, C0):
x = x - t * dx / C0
x = np.mod(x, h + 1e-9) # hack to avoid shift in the plot
return np.array((0.1 < x) * (x < 0.3), dtype = "float")
solve(gauss, 1, 200, "task02_gauss_C0_1", analytic = gauss_analytic)
solve(square, 1, 200, "task02_square_C0_1", analytic = square_analytic)
solve(gauss, 0.7, 200, "task02_gauss_C0_0.7", analytic = gauss_analytic)
solve(square, 0.7, 200, "task02_square_C0_0.7", analytic = square_analytic)
solve(gauss, 1.1, 200, "task02_gauss_C0_1.1", analytic = gauss_analytic)
solve(square, 1.1, 200, "task02_square_C0_1.1", analytic = square_analytic)
# To generate videos out of the generated plots
# $> ffmpeg -y -r 60 -i plots/task02_gauss_C0_1_%05d.png -pix_fmt yuv420p task02_gauss_C0_1.mp4
# $> ffmpeg -y -r 60 -i plots/task02_square_C0_1_%05d.png -pix_fmt yuv420p task02_square_C0_1.mp4
# $> ffmpeg -y -r 60 -i plots/task02_gauss_C0_0.7_%05d.png -pix_fmt yuv420p task02_gauss_C0_0.7.mp4
# $> ffmpeg -y -r 60 -i plots/task02_square_C0_0.7_%05d.png -pix_fmt yuv420p task02_square_C0_0.7.mp4
# $> ffmpeg -y -r 60 -i plots/task02_gauss_C0_1.1_%05d.png -pix_fmt yuv420p task02_gauss_C0_1.1.mp4
# $> ffmpeg -y -r 60 -i plots/task02_square_C0_1.1_%05d.png -pix_fmt yuv420p task02_square_C0_1.1.mp4

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

BIN
Exercise_03/videos/task02_gauss_C0_0.7.mp4 (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/videos/task02_gauss_C0_1.1.mp4 (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/videos/task02_gauss_C0_1.mp4 (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/videos/task02_square_C0_0.7.mp4 (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/videos/task02_square_C0_1.1.mp4 (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/videos/task02_square_C0_1.mp4 (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/videos/video_1_1.mp4 (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/videos/video_1_2.mp4 (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/videos/video_1_3.mp4 (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/videos/video_1_4.mp4 (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/videos/video_1_x.mp4 (Stored with Git LFS) Normal file

Binary file not shown.

BIN
Exercise_03/videos/video_2_gauss_C0_1.mp4 (Stored with Git LFS) Normal file

Binary file not shown.