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
4 个评论
Torsten
2025-7-14
As a sidenote: the equation is called "Telegraph Equation".
William Rose
2025-7-14
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. .
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. . 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
2025-7-15
One can interpret the equation 
as the wave equation for an elastic membrane with fluid resistance and uniform-in-space external forcing, F(t).
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.
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.采纳的回答
William Rose
2025-7-14
编辑：William Rose
2025-7-14
[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.
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.
5 个评论
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;
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
2025-7-15
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.
, 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
2025-7-15
移动：Torsten
2025-7-15
William Rose
2025-7-15
@Muhammad Fahim, you're welcome.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
