Help me in Solving PDE

How to solve the below problem in MATLAB with finite difference method.

equation:

λ ∂²w/∂t² + ∂w/∂t = Re² e^(iωt) + ∂²w/∂x² + ∂²w/∂y²

initial conditions:

w(t=0) = 0, ∂w/∂t(t=0) = 0

boundary conditions:

w(x=±1) = 0, w(y=±r) = 0

parameters:

λ=6, Re=10, ω=1, r=1, t_final=pi/2.

Stability condition: dt<=(dh^2)/8, where dh=dx=dy.

The figure for these values is attached here

[Torsten]

As a sidenote: the equation is called "Telegraph Equation".

[William Rose]

@Torsten,

Thank you for pointing that out. I've encountered the version with one spatial dimension as a model of impulse transmission along nerve axons. I haven't seen it with two spatial dimensions before. In the version posted here, the w(x,y,t) term (with no partials) is absent. This happens in the telegraph equation if it is a lossless transmission line (R=G=0). In that case, the term and the term both vanish. But in the problem here, the term does not vanish. So it's not a totally lossless transmission line (or surface, since it has two spatial dimensions): there are series losses but not parallel losses (R>0, G=0), or there are parallel losses but not series losses (G>0, R=0). At least that's what I make of it. .

[Torsten]

@William Rose

Thank you for the background information. For me as a mathematician, the names of the equations usually help me to find adequate discretization schemes from the literature. Many times I wished I had a better understanding about the physical background - it would have made it much easier to interprete and occasionally discard results.

[William Rose]

@Torsten,

One can interpret the equation as the wave equation for an elastic membrane with fluid resistance and uniform-in-space external forcing, F(t).

The term is the dissipative term due to resistance to motion. If an elastic membrane were actually a fine elastic mesh, then the faster the mesh goes through the air, the more drag there is on it.

Answer 1

[William Rose]

@Muhammad Fahim,

[Edit: fix typos (none in code).]

Write the finite difference equations and solve it explicitly. This approach is not very elegant or efficient, but it will work if you choose and and small enough, and so they satisfy the stability criteria.

Your equation includes Re² e^(iωt) , where Re=10. Are you trying to take the real part, or do you really want a complex driving force and a complex solution?

Let w be indexed by i,j,k, which correspond to x, y, and t respectively. Then the finite difference equation can be written

Rearrange the equation above to solve for w(i,j,k+1), which is w() at the next time point:

Rearrange some more:

The equation immediately above is what you want: an equation for w at the next time point, in terms of w at earlier time points.

Write code with three nested for loops. The outside loop is over k (time). The inside loops are over i and j (x and y). In the inner-most loop, update w. Enforce the temporal and spatial boundary conditions.

Here is some pseudo code:

% define dx, dy, dt appropriately, respecting the stability criteria
% define imax, jmax, kmax appropriately
w=zeros(imax,jmax,kmax) % initialize w
for k=2:kmax-1
    for i=2:imax-1
        for j=2:jmax-1
            w(i,j,k+1)=(right hand side);  % use right hand side of long equation above
        end
    end
end

k=1,2 corresponds to t=-dt, t=0. w=0 at these times, according to the initial conditions, and we have initialized it as such. The code computes the values for w(i,j,k) starting at k=3, which corresponds to t=+dt. This allows me to apply the initial conditions without getting an invalid array index error.

The i loop goes from 2 to imax-1 so that there won't be an array index error at the edges, and this also enforces the boundary condition on x, which is that w=0 at x=+1,-1.

The j loop goes from 2 to jmax-1 so that there won't be an array index error at the edges, and this also enforces the boundary condition on y, which is that w=0 at y=+1,-1.

Check my equations above carefully, to make sure there are no mistakes.

Good luck.

[Muhammad Fahim]

@William Rose

I developed the below FDM scheme, but it did not generates the graph that I have attached here.

Re = 10; lambda = 6; omega = 1; r = 1;
xmin = -1; xmax = 1; ymin = -r; ymax = 1;
dh = 0.02; dt = (dh^2)/8; Tfinal = pi/2;

Nx = round((xmax-xmin)/dh); 
Ny = round((ymax-ymin)/dh);
x = linspace(xmin,xmax,Nx); 
y = linspace(ymin,ymax,Ny);
Nt = round(Tfinal/dt);

w = zeros(Nx,Ny,Nt);
% Initial conditions
w(:,:,1) = 0;  % Initial condition at t=0
w(:,:,2) = 0;  % Initial condition at t=dt

for k = 2:Nt-1
    t = k*dt;
    for i = 2:Nx-1
        for j=2:Ny-1
            dpdz = Re^2*exp(1i*omega*t);
            wxx = (w(i+1,j,k)-2*w(i,j,k)+w(i-1,j,k))/(dh^2);
            wyy = (w(i,j+1,k)-2*w(i,j,k)+w(i,j-1,k))/(dh^2);
            a11 = (lambda/dt^2)+1/(2*dt);
            b11 = (lambda/dt^2) - 1/(2*dt);
            c11 = (2*lambda)/(dt^2);
            w(i,j,k+1)=(c11*w(i,j,k) - b11*w(i,j,k-1) + dpdz + wxx + wyy)/(a11);
        end
    end
    % Boundary conditions
    w([1 Nx],:,k+1) = 0; 
    w(:,[1 Ny],k+1) = 0;
end

[X, Y] = meshgrid(x,y);
figure(1);
surf(X,Y,real(w(:,:,end))');
xlabel('x'); ylabel('y'); zlabel('w');
title('Solution at final time');
shading interp;

●

[Torsten]

@Muhammad Fahim

Check your coefficients a11, b11 and c11 - I think they are not correct.

I get

a11 = lambda/dt^2 + 1/dt;
b11 = lambda/dt^2;
c11 = 2*lambda/dt^2 + 1/dt;

Do you have a reference to a publication / code where the image was presented ?

PS:

I see now that you discretized dw/dt as (w(i,j,k+1)-w(i,j,k-1))/(2*dt) instead of (w(i,j,k+1)-w(i,j,k))/dt - so your code should be ok.

[William Rose]

@Muhammad Fahim,

The surface plot you have presented looks reasonable. I say that because the external forcing function,, which is the same everywhere on this square surface, has , and Tfinal=pi/2. Therefore the external force goes through one quarter of a cycle during the period of integration, and its rate of change is zero at t=Tfinal. The spatial boundary conditions are that the membrane is clamped to zero at the edges x=+-1, y=+-1.

[Muhammad Fahim]

@William Rose @Torsten

Thank you so much for your help. I truly appreciate your time.

[William Rose]

@Muhammad Fahim, you're welcome.