Solution scheme for solution to PDE independent of time step

[robby991]

Hi, I am using the central difference method to solve a diffusion-based partial differential equation. However, my code now will not run because the time step has to be so large that Matlab cannot handle it. The large time step is due to the stability of

stability:

 Diffusion coefficient*time step/(space step^2) < .5

Where my space matrix goes from 0 to X with spacings of dx as determined by the number of space nodes (numx) as

dx = x/numx-1

my time matrix goes from 0 to t with spacings of dt as determined by the number of time nodes (numt) as

dt = t/numt-1

The diffusion coefficient is fixed, my space steps are essentially fixed because I need to have a certain amount of space nodes in my mesh so I cannot reduce them. My total time cannot go below a certain value either for my code to make sense, so the only thing I can do is increase the number of time steps to make dt as small as possible. However this number is unreasonably large that Matlab cannot handle it and runs out of memory.

Is there any other solution method which is essentially independent of time stepping; as there is no stability that requires on the size of the time step. Any computer differential equation solution packages? Thank you.

 

[AlephZero]

"Diffusion" probably means your PDE is parabolic. In that case Crank-Nicholson type methods are unconditionally stable and widely used in practice.

Of course might still be a relation between the spacing between the mesh points and the time steps for good accuracy (i.e. if you have a finer mesh, you need smaller time steps to resolve what is happening), but that comes from the behaviour of the solution of the PDE, not from the limitations of the numerical method.

 

[Chestermiller]

Hi, I am using the central difference method to solve a diffusion-based partial differential equation. However, my code now will not run because the time step has to be so large that Matlab cannot handle it. The large time step is due to the stability of

stability:

 Diffusion coefficient*time step/(space step^2) < .5

Where my space matrix goes from 0 to X with spacings of dx as determined by the number of space nodes (numx) as

dx = x/numx-1

my time matrix goes from 0 to t with spacings of dt as determined by the number of time nodes (numt) as

dt = t/numt-1

The diffusion coefficient is fixed, my space steps are essentially fixed because I need to have a certain amount of space nodes in my mesh so I cannot reduce them. My total time cannot go below a certain value either for my code to make sense, so the only thing I can do is increase the number of time steps to make dt as small as possible. However this number is unreasonably large that Matlab cannot handle it and runs out of memory.

Is there any other solution method which is essentially independent of time stepping; as there is no stability that requires on the size of the time step. Any computer differential equation solution packages? Thank you.

This is a classic example of a set of stiff ODE's, produced by applying the Method of Lines to a parabolic PDE. A standard way of dealing with stiffness of the ODEs is to go implicit, and solve for the dependent variables simultaneously at time t + Δt. To do this you need to switch from the Forward Euler integration formula you have been using to the backward Euler formula:

Forward Euler:
x (t+ Δt) - x (t) = F (x (t )) Δt


Backward Euler:
x (t+ Δt) - x (t) = F (x (t + Δt )) Δt

Using the backward euler formula forces you to solve a set of linear algebraic equations at each time step. But it allows you to take as large a time step as you want, with very good accuracy and no instability.

This method of solving sets of stiff ODEs was perfected upon and expanded upon by C. W. Gear in the 1960's. It is probably one of the biggest advances in numerical analysis ever. Commercial and free ware automatic integration Gear Packages are readily available. They are very powerful. Using these packages frees you from having to do all the details of the calculations. See if MATLAB has a stiff ODE solver option.

 

[robby991]

Thank you both for your response. Based on your feedback, I am doing the crank nicholson implicit technique. I can code the diffusion equation in Matlab with the Crank Nicholson scheme using backslash divide on a martix, but there is an added term to the end of my equation that I do not know how to fit it into the tridiagonal martix. See the attachment, there is a term VF/K+F added to the end of the equation (where we are solving for F). The use of a tridiagonal matrix to solve this equation requires linear operation, that is

AF(t+dt) = BF(t) +d(t)

How do I rearrange that last equation after the diffusion equation (will probably be t, not t+dt?), so the F from the denominator is now in the numerator, so I can incorporate it into the matrix? I am unsure how to solve this, it is really puzzling me. If I can understand how to set up the matrix I can code it no problem. Thanks again.

 

[robby991]

Does anyone have any idea how this V*F/K+F would be "linearized" and fit into the matrix?

 

[bigfooted]

Normally you would take an average of the current and next timestep:
$\frac{VF}{K+F} \approx \frac{1}{2}(\frac{VF^{n+1}}{K+F^{n+1}} + \frac{VF^{n}}{K+F^{n}})$
This poses a problem because it's nonlinear. The simplest approach is to assume it's almost constant:
$\frac{VF}{K+F} \approx \frac{VF^{n}}{K+F^{n}}$
But this is not always stable.
You can also use a two-step method: use this solution as an estimate of $F^{n+1}$ and calculate the nonlinear term on n+1. Then use the first method above with a known nonlinear term.

More advanced is to combine this with a Newton method. For the Jacobian J, you can derive the exact Jacobian of the Crank-Nicolson method or use an approximation.
You can also switch to an explicit method like Runge-Kutta.

 

[robby991]

Thank you for your response. I actually have been trying your second method, where j = time steps:


$\frac{Vmax*F(j+1)}{KM+F(j)}$


However, I guess my question is how to actually solve the matrix. If I set up my matrix of the form

A*F(j+1) = B*F(j)

then,


$F(j+1)*\frac{Vmax}{KM+F(j)}$

ends up appears on the left side of the equation, along with the rest of the F(j+1). However, it contains a F(j) term in the demonimator that becomes part of "A". Therefore, it cannot be solved in Matlab, by say the "backslash divide" technique. Am I correct? Shouldn't all the F(j+1) appear on one side and F(j) sppear on the other? I can solve the diffusion equation using crank nicholson and the backslash divide in matlab, but I am not understanding how this equation now fits into the matrix.

If someone can clearify this that would be much help, or propose a different technique.

 

[bigfooted]

You need to construct the tri-diagonal matrix only for the quantities on the new time step n+1. You calculate the nonlinear term using the known values on the old time step n. The nonlinear term is therefore known and stays on the right-hand side.