We changed our privacy policy. Read more.
10

I have found that I must keep the value of dt near dx or the results become unstable This behavior you have noticed is known as the Courant–Friedrichs–Lewy (CFL) condition. Indeed, there are timestep restrictions that are related to spatial refinement and choice of time integration method. Due to the diffusion term in 13a, the timestep must satisfy an ...


8

Solving partial differential equations with explicit timestepping methods relies on meeting a certain CFL condition for stability. Since you are using a forward Euler timestepping scheme, you must ensure that your timestep is small enough compared to the grid spacing. In the case of the heat equation (Fick's 2nd law is a heat equation), your CFL number $C$ ...


8

It's easy to derive that equation from Fick's law. You have this diffusion equation as: $$\frac{\partial C}{\partial t} = D \nabla^{2} C$$ The mean square displacement weighted by the concentration profile is defined as: $$\langle r^{2}(t) \rangle = \int_{\Omega} |\vec{r}|^{2} C(\vec{r},t) d^{3} \vec{r}$$ The time-evolution of this mean square ...


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

This is possible (see [1]) but uncommon, as it requires Monte Carlo moves that alter the current conformations by a very small perturbation. In that setting of "small" Metropolis MC moves, it is usually easier (both in theory and in practice) to just use Molecular Dynamics instead. [1] Kikuchi, K., Yoshida, M., Maekawa, T., & Watanabe, H. (1991). ...


5

It is in some cases possible to map the dynamics obtained in MC simulations to other (more realistic) dynamics, especially for the case of dense colloidal suspensions. The following two papers talk about the problems and caveats of performing such a mapping: http://aip.scitation.org/doi/10.1063/1.3414827 (spherical particles) http://aip.scitation.org/doi/10....


5

First off, the PDE can be rewritten instead as $$\frac{\partial C}{\partial t} = \frac{\partial}{\partial x}C\frac{\partial C}{\partial x}$$ or, by applying the product rule in reverse again, as $$\frac{\partial C}{\partial t} = \frac{1}{2}\frac{\partial^2}{\partial x^2}C^2.$$ This equation is often referred to as the porous medium equation or the slow ...


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

You are starting from a uniform temperature and you have insulated boundary conditions; so there is no heat conduction occurring. Likewise, your initial mole fraction is also constant in $x$ so that the heat input is uniform along the length. So the fact that the temperature and mole fraction don't change as a function of $x$ is exactly what you should ...


4

That's a modelling problem, not a CS one. This will depend on what is diffusing and on the reasons why 2D is relevant. You may want to give physical details and ask this on physics.SE. Alternatively if you just want to get your hand in with solving diffusion, then the numbers don't have to reflect a physical fact (especially if you have diffusion only, the ...


4

More generally, the equation you are solving is \begin{equation} \frac{\partial u}{\partial t} = \nabla \cdot \left( D(x) \nabla u \right) \quad , \end{equation} where the first nabla operator represents the divergence operator, and the second represents the gradient operator. Your equation is the special case for a one-dimensional cartesian coordinate ...


3

You just add the diffusion along the other dimensions. This superposition from orthogonal directions makes some sense, as they are independent. So, going by wikipedia for Fick's second law of diffusion in 1D: $$ \frac{\partial \psi}{\partial t} = D \frac{\partial^2 \psi}{\partial x^2} $$ We extend it to 2d as: $$ \frac{\partial \psi}{\partial t} = D \...


3

In statistics the widely used test for checking if the distribution is gaussian is the Jarque-Bera test. Koizumi [1] presents an equivalent test for the multivariate case. I don't know if there something ready to use in C++, Fortran, or Python. However, at first sight, the test looks easy to implement. [1] - Koizumi, Kazuyuki, Naoya Okamoto, and Takashi Seo....


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

There is no reason to believe that two random fields with the same arithmetic mean would yield solutions that have anything to do with each other. In fact, for the case you consider, one might imagine that maybe the harmonic mean is actually a better indicator, but even that is unclear -- it could also be the geometric mean. Apart from this, you have to ...


2

You should never use explicit method for the diffusion equation. Implicit is unconditionally stable and just as easy to implement. Also if you use an implicit method (like backward Euler or Crank-Nicolson) it will not matter how small d is. In fact you could use a dirac delta function if you wanted. As far as speed goes, doing a implicit method with a ...


2

For the heat equation, the wave speed is infinite and this is reflected by the fact that you are solving a globally coupled differential equation if you are using an implicit time stepping scheme.


2

You need to unfold the positions after the simulation is done. I use the following subroutine: subroutine unfold_positions(L, X, Xu) ! Unwinds periodic positions. It is using the positions of particles from the ! first time step as a reference and then tracks them as they evolve possibly ! outside of this box. The X and Xu arrays are of the type X(:, j, i), ...


2

Your first test should be to compute the mean value and covariance matrix of your point sample. If these converge to the correct values as you increase the number of samples, you are, from a practical perspective on the safe side that your points indeed come from the correct distribution. Of course, in practice there are many distributions that have the ...


2

I would do it exactly how you describe, with the two following remarks. I would solve the two evolutions simultaneously. In the way you describe, you will need to store and retrieve all the time steps for the solution $u_A$. Instead, you could solve sequentially for $u_A$ and $u_B$ at each time step, since $u_B$ just requires two solutions for $u_A$. If you ...


2

At first glance, it looks like you are using the method of lines with forward Euler time steps. What WolfgangBangerth is getting at with his example is that even for a simple heat equation, the stability limit of forward Euler (namely, that $|\lambda\Delta{t}| < 1$) combined with the eigenvalues induced by a finite difference approximation of the ...


2

Interesting question. I too would like to know if this is mentioned anywhere explicitly. My guess is maybe this is a little too hard and not too useful to be used/taught. Perhaps someone can give a proper reference. My calculations (see code below) suggest, if I am not wrong, that for the equation $$ u_t = \mathcal{L}u, \qquad \mathcal{L}u = \nabla\cdot(D\...


2

One of the best ways to test a PDE solver is to use the method of manufactured solutions. Essentially, you modify the PDE (and discretization) by adding a source term that yields an exact solution known in advance. You can then compare your numerical solution against the exact solution for debugging purposes. You should test your numerical solution against a ...


2

It's a standard advection-diffusion equation. As long as your coefficients are bounded away from zero, there is really no difficulty associated with this equation with the possible exception of the fact that you need to stabilize the advection term if it dominates the diffusion term. I would use a standard finite element method, plus something like SUPG for ...


2

I assume you're trying to solve an equation that looks like: \begin{align} -\nabla \cdot (a(x)\nabla{u}) = f, \end{align} for $x$ in some domain $\Omega$, although the same approach would be fine (for a residual evaluation, anyway) if $a$ were also a function of $u$. The stiffness matrix will take the form \begin{align} A_{ij} = \int_{\Omega}a(x)\nabla\...


2

You are applying Richardson extrapolation to increasingly-accurate independent sample paths of an OU process. How would this work even for independent identically normally distributed random variables $y_1,y_2,y_3,y_4\sim\mathrm{N}(\mu,\sigma^2)$? These have already "converged" to the true distribution $\mathrm{N}(\mu,\sigma^2)$, and you might expect the ...


2

If you want to obtain the diffusion coefficient for 2D from 3D you have to take a look how the 2D diffusion equation is derived from 3D equation. It will be integrated with respect to one of the spatial dimensions. Diffusion coefficient is a $3 \times 3$ tensor in 3D, in 2D it will be $2 \times 2$ tensor. So in my opinion you cannot simply take the value ...


2

In fact, your equation is a non-linear advection-diffusion. Due to the fact that your problem is time-dependent, it could be easily solved by finite-difference: $$\frac{\partial z}{\partial t} = -C (\sin(\omega t))^{m} x^{hm} n (\frac{\partial z}{\partial x})^{n-1} \frac{\partial^{2} z}{\partial x^{2}} -C (\sin(\omega t))^{m} x^{hm} h m x^{-1} (\frac{\...


2

It seems still you don't specify your boundary conditions explicitly despite the suggestion given in one of your previous questions. As far as I understand from your MATLAB code, your boundary conditions don't make sense at all. You have two Neumann boundary conditions at the left and right sides of your 1D domain: $$-D\frac{\partial C}{\partial x}|_{x=0} = ...


2

Short Answer: There is no general result that would hold for all implicit schemes. The reason is that how your method behaves with respect to numerical diffusion depends on the specific combination of numerical time stepping scheme and spatial finite difference. However, note that if your discretisation is consistent with your PDE and your PDE does not have ...


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