People who code: we want your input. Take the Survey

Tag Info

20

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} \... 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. 6 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 Yes, the problem with mixed boundary conditions is well posed. What's not clear to me is this: Why do you approximate the derivative via the two-sides approximation? Shouldn't it be enough to just the following?$$ \frac{f_6 - f_{6-\Delta x}}{\Delta x} = 0. $$In your animations, what is the size of \Delta x? Your curve looks very smooth, which can be ... 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 ... 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 ... 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 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 ... 3 Crank-Nicholson is usually used for linear ODE but can also be extended to a nonlinear ODE$$\dot z = f(z)$$like so:$$\frac{z_{n + 1} - z_n}{\delta t} = \frac{f(z_{n + 1}) + f(z_n)}{2}.$$Another similar alternative is the implicit midpoint rule:$$\frac{z_{n + 1} - z_n}{\delta t} = f\left(\frac{z_{n + 1} + z_n}{2}\right).$$The Crank-Nicholson and ... 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 ... 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? 2 User 03161 asserts that the Crank Nicolson method is not appropriate for advection problems, but boyfarrell provides a working code with results visualized in a movie. In fact they are both correct, but neither gives the full perspective. At the beginning of boyfarrell's answer the correct C-N formula for linear advection is written out. In that formula are ... 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 It becomes a root finding problem. Find u^{n+1} such that$$ u^{n+1} + \Delta t \theta F(u^{n+1},t^{n+1}) = u^n + \Delta t (\theta - 1) F(u^{n},t^{n}).$$Now, you can multiply by a test function, form your trilinear form w(\cdot,\cdot,\cdot) to obtain the weak problem; Find u^{n+1} in an appropriate space V such that$$w(u^{n+1},u^{n+1},v) = (u^n +...

1

Take the first expression and start to reduce the $x$ values to $x_i$, \begin{align} \frac{u_{i+1}^n - u_i^n}{x_{i+\frac{1}{2}}} - \frac{u_i^n - u_{i-1}^n}{x_{i-\frac{1}{2}}} &= \frac{(x_i-\frac12Δx)(u_{i+1}^n - u_i^n) - (x_i+\frac12Δx)(u_i^n - u_{i-1}^n)}{x_{i-\frac{1}{2}}x_{i+\frac{1}{2}}} \\ &=\frac{x_i}{x_i^2-\frac14Δx^2}(u_{i+1}^n - 2u_i^n + u_{...

1

Problem well posed Your problem is well posed. On the discretization with Crank-Nicholson I am not familiar with MMS, and I wonder how you got that ungeneralised form of the diffusion equation. Anyways, as I understand it, you are using the Crank-Nicholson method to discretize the following differential equation: $\frac{\partial f}{\partial t} = x^2 \... 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

I was curious about your problem and it was not difficult for me to test CN method for your two equations. Added - parts beginning with "added" were written afterwards to address the comment. The point given by David is important. Here are three plots for initial condition (added - that do not fulfill the suggested condition $\alpha u_{xx} + \beta v_{xx}=0$)...

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 ...

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