8

One situation in which the usual method-of-lines approach cannot be used in a straightforward way is with equations that have mixed space-time derivatives.. By "usual method-of-lines approach", I mean discretization of spatial derivatives followed by application of a Runge-Kutta or linear multistep method. This usually applies only to systems of first-...


7

To answer your questions in order: Any implicit method for solving an ordinary differential equation involves solving a system of nonlinear equations. You can do this through variants of Newton's method, successive substitution, full approximation scheme, or any other approach that solves systems of nonlinear equations. (The caveat is, of course, depending ...


6

You are correct: If you satisfy the CFL condition, then all that guarantees is that your scheme is stable, i.e., the numerical solution does not go to infinity. But the CFL condition says nothing about how accurate the numerical solution is. For that, indeed $\Delta z$ and $\Delta t$ must also be small enough compared to the features of the exact solution. ...


5

I have two extra points I would like to add to Wolfgang's answer. A formulation of the CFL condition that I find more useful than the classic formula is this: A necessary condition for the stability of a numerical scheme is that the numerical domain of dependence bounds the physical domain of dependence. This is exactly what good old $$ \dfrac{\Delta ...


5

A nice reference for this is Chapter 10 of LeVeque's finite difference book. Of course, it only covers basic finite difference approaches, and there are plenty of others (all within the method of lines framework). The method of lines is indeed applicable in multiple dimensions, and two-dimensional problems are discussed in the reference just given. Most ...


5

A timestepping method can be considered always as an iterative solver for the steady state problem if the solution is guaranteed to run into a steady state. There are several aspects here: If the timestepping scheme is implicit and sufficiently stable, you can increase the time steps when you approach the limit. This is often referred to as pseudo ...


5

ode15s is designed to handle stiff systems of ODEs so I doubt if the problem you are encountering is that your "equations are too stiff" It is more likely that your spatial discretization has an error for some reason or your have some other MATLAB programming error. I suggest the following approach to debug this: Set the final integration time for ode15s to ...


4

There are many definitions of DAE index. MATLAB is probably referring to the differentiation index. A few integrators can solve Hessenberg form index-2 DAEs (I believe IDA can do this, along with a few other packages), but most require the index to be 1 or less. Reducing the index of your DAE requires manipulating the underlying equations in the DAE, e.g., ...


4

It is easier to see what is happening if you put the $\partial B/\partial z$ term also on the right hand side: Then, you have one equation with a time derivative and one that does not. This is what we typically call a "Differential Algebraic Equation" (DAE), and it is the same as what is happening, for example, in the time dependent Stokes and Navier-Stokes ...


3

If your advection problem had Dirichlet of Neumann boundary conditions, the linear system would be tridiagonal and you could apply the Thomas algorithm. With periodic boundary conditions, however, we lose this. If c(x) is a constant independent of x, the matrix would be circulant and linear systems could be solved efficiently using FFTs. An even better ...


3

All of the MATLAB ode solvers adjust the integration step size to try to keep the error in the solution less than certain prescribed tolerances. If the error exceeds these tolerances, the step size is reduced until it reaches essentially zero (1.734723e-18). Usually this error message arises because the problem is ill-posed. If you terminate the solution at ...


2

This is a fairly big question. The basic idea is that in solving $\partial_t u = \mathcal{L}u$, you approximate $\mathcal{L}$ with a matrix $A$, and solve $$ u' = Au. $$ In this formulation you have two main issues: What method to use to solve the ODE. The method has to be implicit because of the nature of parabolic PDEs, and the stiff systems of ODEs that ...


2

If your discretization of the reaction part of the problem is explicit, you can still use Prof. Bangerth's suggestion. You can quickly compute $K*u$ using the fast Fourier transform on both $K$ and $u$, multiplying the transforms and then transforming back. Then you can easily compute $u\cdot K*u$ at any grid point. However, you may run into stability issues ...


2

Do a Fourier transform in space and solve the ODEs that describe the first $N$ Fourier coefficient. This is possible because the Laplacian turns into a multiplication with the wave number upon Fourier transform, and the convolution turns into a straight multiplication. This makes all operations reasonably trivial, particularly if all you have is a 1d domain.


2

What you describe is called Differential-Algebraic Equation (DAE) system. Depending on the index of the system, these can be easy to solve or very hard. Take a look at: http://www.mathworks.com/matlabcentral/fileexchange/7481-manuscript-of--solving-index-1-daes-in-matlab-and-simulink- Solvers ode15s and ode23t of Matlab can handle index-1 DAEs. If the ...


2

Let's take the example of the unsteady one-dimensional heat equation inside a solid on a domain $x\in[0,1]$: $$\partial_t u - D\partial_{xx} u = 0$$ with the initial profil $u(0,x) = u_0(x)$ at $t=0$, and Dirchlet boundary conditions enforcing that the wall at $x=0$ (respectively $x=1$) is at temperature $u_{L}$ (respectively $u_{R}$). If we use discretise ...


2

You can convert your equation into the following system of two PDE $$ \frac{\partial u}{\partial t} = v $$ $$ 0= \frac{\partial^2 u}{\partial x^2} + A \frac{\partial^2 v}{\partial x^2} - B $$ There are now two dependent variables, $u$ and $v$. The discretization in the spatial dimension (i.e. method of lines (MOL)) is straightforward; central difference ...


2

Your question is not quite clear, but if your problem is roughly analogous to $\frac{d\mathbf{x}}{dt}=\mathbf{f}[\mathbf{x},y]$ , $g[\mathbf{x},y]=0$ then you just need to use a mass matrix in your call to ode15s. This is how you "tell" Matlab that you have an algebraic equation. That is, in Matlab pseudo-code you can express the above as eye(n)*dxdt = f(...


2

What do you mean by "I have realised that several things did not affect the mass conservation". How did you measure conservation ? Are you asking if your method conserves total mass ? My answer below is based on a conservation law I notice in your model. You wrote an $r$ in your PDE but did not use it later, so I am assuming it is zero. Your model has this ...


1

It is useful, as a first step, to analyze method-of-lines (MOL) treating the time integration as exact. For example, suppose we are solving the 1D advection equation $ \frac{\partial{n}}{\partial{t}} =-c \frac{\partial{n}}{\partial{x}}, $ where $c$ is the advection speed. Then, for example, for central difference, we'd have $ \frac{d}{dt} n_i = - \frac{c}{2 ...


1

I think you are using a downwind- instead of an upwind finite difference. This leads to your code imposing a boundary condition where it is not allowed. The solution to your convection equation is basically (ignoring the left BC for the moment) $$ C(x,t) = C_0(x - v t) $$ where $C_0$ is your initial value. Thus, if $v > 0$, it is a rightward travelling ...


1

It is easy to see that the function \begin{equation} P(t) := \int \limits_0^1 p(x,t) \, \mathrm{d}x, \quad t \geq 0, \end{equation} is constant over time (i. e. $P(t) = P(0)$, $\forall \, t \geq 0$), if $p$ is a solution of your initial-boundary-value problem. If you discretize the problem in space as you described, then you obtain a system of linear ...


1

If there are only continuous solutions in your problem, the finite difference method is enough. However, if there exists shock, the finite difference/finite volume WENO/ENO suit for your purpose. Generally speaking, you should first discretize the systems in space and make it to a systems of ODEs. Then, use Runge-Kutta methods to evolve in time.


1

Please try to improve your explanation: with the index i do you refer to the space variable or time one? Anyway, as far as I understood from your post, I think that you need to enforce some boundary condition if you're dealing at the end/beginning of your spatial domain. Usually, according to the physics of your problem, you have to enforce Neumann or ...


1

In short, yes, you can solve the entire semi-discretized PDE at once using ode15s, or any other ODE solver (ignoring, for the moment, issues of stability and accuracy). I would not arrange the unknowns from the semi-discretized system as a matrix; instead, I would arrange them as one long vector. Although indexing this vector will be cumbersome, expressing ...


1

You could use a stencil to do a multidimensional least-squares fit which will account for transverse components. I've seen this done in finite volume before, [Weller et al. 2009, Weller & Shahrokhi 2014] give details of using upwind-biased stencil on unstructured two-dimensional meshes, but the same techniques can be straightforwardly applied to uniform ...


1

Here's a 30 minute introduction that explains the difference between the method of lines and the Rothe method: http://www.math.tamu.edu/~bangerth/videos.676.26.html


1

For reference, if you want to learn about the difference between the Rothe method and the method of lines, then maybe lectures 26 and 27 on http://www.math.tamu.edu/~bangerth/videos.html are of interest to you. If you want to see how to actually implement them, take a look at lecture 29. (Disclaimer: these are my own lectures.)


1

I don't have enough reputation to comment, but I hope this answer helps you. Personally, I find that I learn best through example. Try looking at the code here to see how MOL was implemented in Python with centered finite difference approximation (an ODE solver was used). If you'd like to learn more about solving ODEs, parabolic PDEs, and MOL, I suggest ...


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