Podcast #128: We chat with Kent C Dodds about why he loves React and discuss what life was like in the dark days before Git. Listen now.
19

Starting with the advection equation is conservative form, $$ \frac{\partial u}{\partial t} = -\frac{\partial (\boldsymbol{v} u)}{\partial x} + s(x,t) $$ The Crank-Nicolson method consists of a time averaged centered difference. $$\frac{u_{j}^{n+1} - u_{j}^{n}}{\Delta t} = -\boldsymbol{v} \left[ \frac{1-\beta}{2\Delta x} \left( u_{j+1}^{n} - u_{j-1}^{n} \...


15

I think that one of your problems is that (as you observed in your comments) Neumann conditions are not the conditions you are looking for, in the sense that they do not imply the conservation of your quantity. To find the correct condition, rewrite your PDE as $$ \frac{\partial \phi}{\partial t} = \frac{\partial}{\partial x}\left( D\frac{\partial \phi}{\...


8

A simplification - the Crank-Nicolson method uses the average of the forward and backward Euler methods. The backward Euler method is implicit, so Crank-Nicolson, having this as one of its components, is also implicit. More accurately, this method is implicit because $u^{n+1}_i$ depends on $F^{n+1}_i$, not just $F^{n}_i$ This is means that the state at ...


8

The maximum principle for Crank-Nicolson will hold if $$\mu \doteq \frac{k}{h^2} \leq 1$$ for timestep $k$ and grid spacing $h$. In general, we can consider a $\theta$-scheme of the form $$u^{n+1} = u^n + \frac{\mu}{2}\left( (1-\theta)Au^n + \theta Au^{n+1}\right)$$ where $A$ is the standard Laplacian matrix and $0 \leq \theta \leq 1$. If $\mu(1-2\theta) \...


6

The algorithm you have implemented is explicit. Crank-Nicolson is an implicit method, and thus requires a solve.


5

I think it might be ill-posed, since the time-dependent parts are linearly dependent. If you add your two time-dependent equations together, you get a time-independent equation: $(\alpha(x)u_x)_x + (\beta(x) v_x)_x = 0.$ Do your initial conditions satisfy this equation?


5

Your equation can be written in the following fashion (any spatial derivative approximation is valid), once space is discretised: $$\frac{1}{c}\frac{du_i}{dt}=-\left(\frac{\partial u}{\partial x}\right)_i(t) + v_i(t) \tag{*}$$ Keep in mind that $v_i(t) = v(x_i,t)$. Now the system of equations depends only on time $t$ you can apply Crack Nicholson method to ...


5

There is something very basic that you should know about hyperbolic problems. Consider the most basic example $\partial_tu+a\partial_xu=0$ with a numerical marching scheme of the form $$u_j^{n+1}=\sum_kc_ku_{j+k}^n$$. This covers all explicit schemes, and all implicit schemes like Crank Nicolson also, if you begin by solving the tridiagonal system. It ...


5

This is exactly the case when the lack of information in the question allows to answer it pretty certainly: it is certainly possible. The error would depend on many factors, including the conditioning of the original problem, particular details of the numerical implementation, and chosen simulation parameters. I do not see any contradiction yet. However, I ...


4

Boundary conditions are often an annoyance, and can frequently take up a surprisingly large percentage of a numerical code. How to implement them depends on your choice of numerical method. Finite difference schemes often find Dirichlet conditions more natural than Neumann ones, whereas the opposite is often true for finite element and finite methods applied ...


4

The solution you believe to be inaccurate is actually by far the more accurate one; you've simply plotted it in a very deceptive way. For $\nu=2$, the exact solution is actually no bigger than about $10^{-35}$ everywhere -- it's zero for all intents and purposes. Therefore the numerical solution is correct to 10 digits -- far better than the accuracy of ...


3

The Crank-Nicolson method is: $\frac{u^{n+1}_{i}-u^{n}_{i}}{dt} = \frac{1}{2}(F^{n+1}_{i}+F^{n}_{i})$ This method calculates the next state of the system, i.e. $u^{n+1}_{i}$, by solving an equation involving the previous states and the next state. In the case of the heat equation for example we would get a linear system and if we are using finite elements ...


3

You should rearrange the terms so that all of the $n+1$ terms are together on one side of the equals sign and all of the $n$ terms are on them other. Then you will have a system of non-linear equations like: $$ A(u^{n+1}_i) u^{n+1}_i=b $$ Where $A$ is a non-linear matrix. Then you can use something like Newton-Raphson to linearize the system and your ...


3

Here is a tutorial on how to solve this equation in 1D with example code. The code is Python (which is similar to MATLAB so you should be able to translate). To extend this to 2D you just follow the same procedure for the other dimension and extend the matrix equation. How to discretize the advection equation using the Crank-Nicolson method?


3

You should reformulate your problem. Let's define the vector $u$ as $u=\left(\begin{array}{c}A\\ B\end{array}\right)$ Then you can write your coupled system as $$\frac{\partial u}{\partial t}=\left(\begin{array}{c}a_0\\ b_0\end{array}\right)\left(u^T \left(\begin{array}{cc}0&0.5\\0.5&0\end{array}\right) u\right)$$ Now we can apply Crank Nicolson as ...


2

If I understand correctly, you are using a centered finite difference in space and the implicit trapezoidal method in time. That scheme is unconditionally absolutely stable, but will generate spurious oscillations. So you should expect to see some increase in the maximum value of $u$, but it shouldn't blow up. If it blows up, you have an implementation ...


1

You need at least one initial condition for $\xi$ and two BC for $\rho$. Where are they? Your equation looks like the heat equation in cylindrical coordinates assuming angular and plane symmetry, with nonlinear heat sources: $$ 2ik\partial_tA-\triangle A=N(A)$$. I would discretise the spatial derivative with finite differences (which is essentialy the same ...


1

$L_x$ and $L_{xx}$ are shorthands (operators) to denote the more extended notation: $$ L_x u_i=(u_{i+1}-u_{i-1}) $$ and $$ L_{xx}u_i = (u_{i-1}-2u_i+u_{i+1})$$. Therefore $L_x$ can be written in the node $i$ as the vector $L_x=(-1,0,1)$ (see the coefficients of $u_{i-1}$ which is $-1$, $u_i$ zero and $u_{i+1}$ one). The same results for $L_{xx}$. For a ...


1

Von Neumann (Fourier modes) stability analysis gives you only a sufficient condition for stability if you compare the amplifying coefficient $r$ with 1. If you have amplifying coefficient bounded by $1 + C\tau$, then after making $\frac{T}{\tau}$ time steps your error will be bounded by $(1 + C \tau)^{\frac{T}{\tau}}$ which has a bounded limit when $\tau \...


1

The scheme is indeed unstable. It explodes - but very very slowly. By printing the maximum eigenvalue of the operator i confirmed the instability. It's greater than 1. Then why does it work? because it's $1.000053263$ and my t_final is small.


1

You seem to have given the 1D equations for the discretizations, even though the problem is in 2D. Regardless, the explicit method requires the least memory since you don't even have to form a matrix to compute the solution at the next time step. If you have the solution vector at time $t_i$, you can apply simple stencil operations and directly obtain the ...


1

it is ok for the discretization. CD scheme has some stability problem when Pe>2, but we can decrease the mesh spacing to obtain a low mesh Pe number. QUICK-scheme is more stable and accurate than CD, so it is ok to implement it. you should increase the order of both convection term and diffusion term to obtain a high order solution. the mesh spacings should ...


1

This question is confusing. At first you are speaking of a steady-state equation, and suddenly you speak of a time scale... I will try to clarify the following. From the numerical PDEs standpoint, you could classify physical phenomena as time-dependent or time-independent. LeVeque's book is a solid reference for this. Usually, time-independent phenomena are ...


1

Let me give an answer that is a general comment on prescribed zero flux for advection-diffusion (or convection-diffusion) PDE that is an important topic and it might be (but not necessary) the problem in your situation. Zero flux boundary condition is very non-standard in this case, so you might check if there is no misinterpretation of your problem. Your ...


1

Thank you for your clarification. The Thomas algorithm solves $Ax=b$ when, e.g., A is a tridiagonal matrix (there are other special cases I believe, but this is not one of them). Your "final form of equation" does not appear to be in this form, nor does it look like this is possible. As @David Ketcheson mentioned, implicit time marching is not very ...


Only top voted, non community-wiki answers of a minimum length are eligible