PDE Toolkit - what is d when m is non-zero?

To solve the following equation

we need to input values of m,d,c,a and f.

I have been able to input m,c,a,f, however, for d, there is a special rule. The following is from Matlab's own support pages.

It says "Generally, d is either proportional to results results.M, or is a linear combination of results.M and results.K". M and K are matrices obtained after the assembleFEMatrices (i.e., after discretization).

Now, to say the least, this is confusing, as I do have the values of the d matrix (2x2 symmetric matrix in my case). Any insights would be appreciated.

[Torsten]

If you have a valid MATLAB licence, I'd contact MATLAB support service.

[Sattik Basu]

Thank you Torsten, I think that is indeed the final option.

Answer 1

[Ravi Kumar]

When m is non-zero in a structural problem, d-matrix (coefficient) could represent a damping matrix. That is the most common use case. Hence, the documentation describes how to compute d, proportional damping matrix, as a combination of global mass and stiffness matrix.

[Torsten]

So independently prescribing d is not possible ?

I think this was the OPs question.

[Sattik Basu]

Ravi, I thank you for taking the time to answer, however, I am afraid the question was misunderstood. 'd' is not necessarily always dependent on 'm' and 'k'. So I wanted to know how one could independently assign those values prior to discretization.

[Ravi Kumar]

Hi Sattik and Torsten,

Yes, specifying 'm' and 'd' simultaneously as coefficients are not supported.

Both 'm' and 'd' typically represent mass matrix, for first- and second-order in time PDEs, respectively. If you are solving a first-order in time PDE then you use 'd', if you are solving a second-order in time PDEs then you use 'm' coefficient. In the latter case, second-order in time, often found in structural dynamics you can specify 'd' as proportional damping matrix of global size.

I understand this is not addressing the issue of specifying both 'm' and 'd' as coefficients and letting the FEM do the discretization. This is a valid enhancement if you can share some details on what does your 'd' represents. Feel free to contact support and provide details.

Regards,

Ravi

[Sattik Basu]

Hi Ravi.

I think it will help if you look at a system where

m = [1 2;

3 4]

d = [5 8;

-8 5];

This is just an example problem. You can consider a,c and f to be arbitrary as well.

The PDE toolkit currently allows me to input m, however, d must be a linear combination. Since I already have a 'd' matrix, help is really appreciated. Thanks

[Torsten]

Maybe setting u3 = du1/dt, u4 = du2/dt and writing the equations as a system of first-order PDEs in time as

du1/dt = u3
du2/dt = u4 
m11*du3/dt + m12*du4/dt + d11*u3 + d12*u4 - terms involving c,a = f1
m21*du3/dt + m22*du4/dt + d21*u3 + d22*u4 - terms involving c,a = f2

is a possibility to circumvent the problem.

[Sattik Basu]

Hi Torsten, this would have worked if the 'au' terms did not exist. The reduction of order would introduce integrals.

[Torsten]

Sorry, but I don't see this.

If you had the PDE for one function u

m*d^2u/dt^2 + d*du/dt - d/dx(c*du/dx) + a*u = f

why not writing it as

du1/dt = u2
m*du2/dt + d*u2 - d/dx(c*du1/dx) + a*u1 = f

or in MATLAB format

[1,0;0,m]*[du1/dt;du2/dt] - d/dx([0,0;c,0]*[du1/dx;du2/dx]) + [0,-1;a,d]*[u1 ;u2] = [0;f]

?

Of yourse I don't know if there are specialized discretization schemes depending on the coefficients prescribed. E.g. if you have the wave equation (m ~= 0), it should be discretized as a hyperbolic PDE, if you have a convection-diffusion equation (m=0), it should be discretized as a parabolic equation.

These specialized schemes could get lost by the first-order transformation.

But all these are questions for the support or the developers of the PDE Toolbox in my opinion.

[Sattik Basu]

@Torsten, your idea about removing m got me into thinking, and if you are willing to write it as an answer to the question I posed, I will accept that answer. I am sure it will help others who are facing similar issues. I was working on the idea you posed, and i was able to build my syste. I have other issues with boundary conditions, but those are future issues.

[Torsten]

Thank you for your feedback.

To see whether the reduction to a first-order system works, I'd test it for the acoustic wave equation

d^2u/dt^2 = omega^2 * d^2u/dx^2

and compare with the solution of the formulation as a PDE second order in time.

[Sattik Basu]

Yes, I am working on validating the answers as we speak. It's always a "state-space" like approach that matlab is comfortable with.