Partial Differential Equation Simulation is giving a solution that goes to infinity

Question

0 votes

Hello, I am trying to simulate a 2D Wave PDE with Dirichlet boundary conditions on 3 of the 4 edges. My issue is that when I try to calculate the evolution of this PDE, the solution reaches values of infinite amplitude, which should not be happening as the motion should be pretty contained. I've done some debugging and i think the issue likely lies in the ODE function (ode15s) used for the method of finite differences. Can someone maybe figure out what is wrong with the code?

%function specification here
function [T,w] = wave2d_test(tf,Nx,Ny,c)
%problem parameters
%c: wave speed parameter
%Nxs: number of nodes in x direction for wave
%Nys: number of nodes in y direction for wave
%Note: Nxs = Nxf & Nys = Nyf
Lx = 1; %length in x direction
Ly = 1; %length in y direction
N = (Nx-2)*(Ny-1); %total number of nodes where w is unknown
%timestepping parameters
t0 = 0; %initial time
%tf: final time
ns = ceil(25*tf)+1; %number of time steps
t = linspace(t0, tf, ns); %time vector
%initialize spatial domain
x = linspace(0,Lx,Nx);
y = linspace(0,Ly,Ny);
[X,Y] = meshgrid(x,y);
dx = x(2)-x(1);
dy = y(2)-y(1);
%initialize solution array W = [wn; vn];
%get initial diplacement and velocity - form W0
wn = winit(X,Y);
vn = vinit(X,Y);
W0 = [wn; vn]; %initial data vector is of length 2*N
%call the timestepping code
tic;
[T,w] = ode15s(@(t,W) RHS(t,W,c,Nx,Ny,N,x,y,dx,dy),t,W0);
odetime = toc;
fprintf('Total run time for ODE integrator: %g\n',odetime);
%postprocess the results
%plot the surface
wmat = [zeros(1,Nx);[zeros(Ny-1,1) reshape(w(end,1:N),[Nx-2 Ny-1])' zeros(Ny-1,1)]];
figure; surf(X,Y,wmat); %plot the final surface
zlabel({'w(x,y)'});
ylabel({'y'});
xlabel({'x'});
title('Solution at final time');
end
function w = winit(X,Y)
Lx = X(1,end); %get Lx
Ly = Y(end,1); %get Ly
[Ny,Nx] = size(X); %get Nx and Ny
%wmat = X.*(Lx-X).*cos(pi*Y/Ly); %define the initial value of w
%wmat = X.*(Lx-X).*cos(pi*Y/Ly).*sin(2*pi*X/Lx); %define the initial value of w
wmat = 0.0*X; %define the initial value of w
figure; surf(X,Y,wmat); %plot the initial surface
zlabel({'w(x,y)'});
ylabel({'y'});
xlabel({'x'});
title('Initial Displacement');
wmat = wmat(2:Ny,2:Nx-1); %trim off left and right and top values
w = reshape(wmat',[(Nx-2)*(Ny-1),1]); %reshape matrix into the w vector
end
function v = vinit(X,Y)
Lx = X(1,end); %get Lx
Ly = Y(end,1); %get Ly
[Ny,Nx] = size(X); %get Nx and Ny
%vmat = (2.0*Y/Ly).*(X/Lx); %define the initial value of w
%vmat = 1.0*sin(pi*X/Lx); %define the initial value of w
%vmat = 0.0*X; %define the initial value of w
vmat = 1.0*Y.*(Ly-Y).*X.*(Lx-X); %define the initial value of w
figure; surf(X,Y,vmat); %plot the initial surface
% Create zlabel
zlabel({'v(x,y)'});
% Create ylabel
ylabel({'y'});
% Create xlabel
xlabel({'x'});
title('Initial Velocity');
vmat = vmat(2:Ny,2:Nx-1); %trim off left and right and top values
v = reshape(vmat',[(Nx-2)*(Ny-1),1]); %reshape matrix into the w vector
end
%define the forcing function f(x,y,t)
function z = f(x,y,t)
z=0;
end
%define the ODE RHS function
function dydt=RHS(t,W,c,Nx,Ny,N,x,y,dx,dy)
fprintf('Time: %g\n',t);
%reshape W into a matrix with 2*(Ny-1) rows and Nx-2 columns
wmat = reshape(W,[Nx-2 2*(Ny-1)])';
w = wmat(1:Ny-1,:); %displacements at interior nodes
v = wmat(Ny:2*(Ny-1),:); %velocities at interior nodes
%initialize dydt
dydt = zeros(size([w;v]));
dydt(1:Ny-1,1:Nx-2) = v; %place v values in the w_t locations
%build the array with ghost nodes
%pad w vector with zeros on all sides
w = [zeros(1,size(w,2)+2); [zeros(size(w,1),1), w, zeros(size(w,1),1)]; zeros(1,size(w,2)+2)];
%compute the RHS of the velocity equation
for j = 2:Ny-1 %which row of nodes
    for i=2:Nx-1 %which column of nodes
        dydt(Ny+j-1,i-1) = f(x(i),y(j-1),t)...
            +(c^2/dx^2)*(w(j,i+1)-2*w(j,i)+w(j,i-1))... %c^2 w_xx term
            +(c^2/dy^2)*(w(j+1,i)-2*w(j,i)+w(j-1,i)); %c^2 w_yy term
    end
end
dydt = reshape(dydt',[2*N 1]);
end