If I understand your PDE it is a linear convection-diffusion-production equation. Convection appear along x1 at a velocity a, diffusion only in x2 with a coefficient b and the production term is linear.

Introducing a spectral representation in space you can solve in time for the Fourier coefficients.

Dear Dr. Filippo Maria Denaro

That was an example, in fact my main goal is to find the stability of Fokker-Planck Equation( convection and diffusion both might appear along x1 or x2), that is a linear parabolic PDE in general, what I really care about is not just the numerical scheme to solve this PDE but also I want a criteria for stability of this kinds of PDE (not numerical scheme stability, but theoretical stability , i.e convergence to a steady state solution regardless of which numerical scheme we use) I find that spectral methods are the basis for finding the stability of the PDEs,

I think I had to solve an eigenvalue problem.

If I understood well, you mean that I had to expand the PDE based on sine or cosine basis functions in general. is it true?

Have a look to the normal mode analysis illustrated in the Chap.12 of the book of Kundu.

Dear Dr. Filippo Maria Denaro

regardless of what method we use (FD or FV) and the order of differentiation approximation (first order, second order, central,...)

finally after discritization we have a linear system

A*U_n+1=U_n

or we can write it

U_n+1=(A)^-1*U_n

were the elements of A matrix is fixed and depends on descritization scheme of derivatives, form of the PDE and its coefficients, source terms and also dt (sampling time).

If we look at the system, it seems it is a discrete dynamical system, that we can talk about its stability based on the eigenvalues of the (A)^-1 matrix.

what do you think about it, do you think we can establish the stability results (theoretical stability not numerical stability) for the PDE based on the eigenvalues of (A)^-1?

First, the general form for a two-step implicit method is

A.un+1 = B.un ->un+1 =(A-1. B).un

The numerical stability can be analyed by the spectral radius of (A-1. B). However, this is a numerical stability analysis, not related to the fluid dynamic stability of the problem.

Dear Dr. Filippo Maria Denaro

you mean that by analyzing the spectral radius of (A-1.B) you can not prove that when t (time) goes to infinity you will reach to a steady state solution?

The matrix governs the evolution of the solution. You can easily see that

u1 =(A-1. B).u0

u2 =(A-1. B).u1 =(A-1. B)2.u0

un =(A-1. B)n.u0

The spectral radius

Dear Dr. Filippo Maria Denaro

I got confused.

I think here I can not understand the difference between the numerical stability and theoretical stability, could you please help me with that.

first consider we have a simple linear ODEsuch as

dx/dt=A*x (1)

with calculating the eigenvalues of A we can prove the stability of the ODE (all eigenvalues are negative)

I am sure that it is the exact definition of theoretical stability i.e existence of a steady state solution.

now we descritize the ODE of system (1) with the first order Euler method

x(k+1)=A*x(k)*dt+x(k) or x(k+1)=(A*dt+I)*x(k),

here if the eigenvalues of (A*dt+I) are in unite circle, ( The spectral radius

Numerical stability and physical stability are two different things.

Talking about the numerical stability there is no reference to a physically steady state but we search for the condition in which the numerical solution is not amplified after repeated application of the discrete time-space (for PDE) scheme.

Classical analysis in matrix form (Gershgorin) or in spectral form (von Neumann) are used to define the conditions (if any) of the numerical stability. For linear scheme, consistence + stability allow us to guaranteee convergence of the solution ( Lax-Richtmyer Equivalence Theorem). You can read the details in textbooks of numerical solutions of PDE.

The physical stability analysis aims to determine the condition in which a physically assessed solution returns to its state after some perturbation is added.

Thank you very much Dr. Filippo Maria Denaro

very nice clarification of the subject.

but intuitively, I think for my example the numerical stability implies convergent to a steady state solution.

I don't see that in your example. If A is a time-independent term, the solution of your ODE is like x(t)=eAt+cost. Steady state, conversely would imply dx/dt->0, that is A->0. In no way I understand the link between numerical and physical stability.

Dear Dr. Filippo Maria Denaro

we do not need A->0

if e^(At)->0 we will also have steady state solution.

and e^(At)->0 when A has negative eigenvalues

Yes, right, I assumed only the positive sign. But as you can see, the condition for the numerical stability implies the constraint on the spectral radius of the iteration matrix, thus no sign is considered.

That is exactly what I seek, some sort of relation between the region of stability of continuous system and the stability region of discrete one, here we can say the half left plane is the stability region for continuous system which mapped into a unit circle for discrete system.

The problem is that the iteration matrix of a discrete system does not necessarily coincide with that of the differential equation.

I mean, while you analyse in a unique way the governing continuous system, you can get different discrete systems simply by changing the discretization scheme. And, consequentely, you can get several cases: conditional stability, unconditional stability or unconditional instability.

Again if we talk about my example, you mean different numerical schemes gives different stability regions?

yes, for example the parabolic heat equation discretized with first order explicit Euler is conditionally stable while the first order implicit Euler is unconditionally stable. That happens using the second order central discretization in space. Changing this latter further changes the numerical stability region.

Dear Dr. Filippo Maria Denaro

I reviewed your answers again and one idea come to my mind that I want to ask you about it.

as you remember we separate numerical stability and physical stability.

In one of your answers you mentioned that :

The physical stability analysis aims to determine the condition in which a physically assessed solution returns to its state after some perturbation is added.

also you mention that :

For linear scheme, consistence + stability allow us to guaranteee convergence of the solution ( Lax-Richtmyer Equivalence Theorem).

what comes into my mind is that if we prove the convergence of our numerical method to the exact solution by (Lax-Richtmyer Equivalence Theorem ), we could be sure that our discrete version of PDE posses the property of the continuous PDE(for example physical stability), so the question is

since our discretization scheme is convergent, can we prove the physical stability of the PDE from its discrete version?

what is your idea about this statement?

could you also please take a look at my other question and its following answers which is related to the above question.

https://www.researchgate.net/post/Is_it_possible_to_prove_the_stability_of_a_continuous_time_system_from_its_corresponding_discrete_version

Unfortunately the things are much more complex. The Lax theorem is valid for linear schemes but if you look to a proof of convergence for the numerical solution of non-linear equations this theorem cannot guarantee convergence.

For a non-linear formulation, Lax and Wendroff theorem states that if a conservative scheme converges it does towards a weak solution. That means we must previously ensure that convergence exists. (see the texbook of Leveque, page 239)

Dear Filippo Maria Denaro

you are correct about nonlinear equations,

but if we choose FD scheme and consider the stability of a linear Fokker Planck equation, is my statement true?

i.e by investigating physical stability of the solution of discretized PDE (where our discretized scheme is converegent based on Lax Equivalence Theorem ) , we could be sure that our Fokker Planck Equation was also physically stable.

Thanks Dr. Denys Dutykh

You mentioned very important fact,

I have some questions and I appreciate if you answer them.

1.

I am not that much familiar with the functional analysis.

for solving the Fokker Planck equation we mainly consider the solution on a bounded domain where on the boundary of the domain we consider the solution to be zero (is this what you mean by the geometry of the domain?)

2. should we choose the basis or we should in fact find the basis by solving the eigenvalue problem for the Fokker planck equation?

3. what do you mean by the completeness of the basis functions?