Help me in Solving PDE

Question

Muhammad Fahim 2025-7-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2178461-help-me-in-solving-pde

评论： William Rose 2025-7-15

采纳的回答： William Rose

picture 1.JPG

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 个评论
显示 2更早的评论隐藏 2更早的评论

Torsten 2025-7-14

编辑：Torsten 2025-7-14

@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 2025-7-15

@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 2025-7-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2178461-help-me-in-solving-pde#answer_1567946

编辑：William Rose 2025-7-14

在 MATLAB Online 中打开

@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

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 个评论
显示 3更早的评论隐藏 3更早的评论

Muhammad Fahim 2025-7-15

移动：Torsten 2025-7-15

在 MATLAB Online 中打开

@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

% 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;

●

Muhammad Fahim 2025-7-15

移动：Torsten 2025-7-15

@William Rose @Torsten

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

William Rose 2025-7-15

@Muhammad Fahim, you're welcome.

请先登录，再进行评论。

Help me in Solving PDE

4 个评论
显示 2更早的评论隐藏 2更早的评论

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Help me in Solving PDE

4 个评论 显示 2更早的评论隐藏 2更早的评论

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论