How can I get Runge-Kutta method to display exponential decay when solving a two dimensional wave equation?
Hello, I am trying to create a wave equation solver using the method of lines and Runge-Kutta to numerically approximate the solution. If my initial conditions are u(0,x,y)=sin(pix)sin(piy), u_t(0,x,y)=0, then I know the exact solution to be e^(-2pi^2)sin(pix)sin(piy). I plot both my approximation (figure 1) and the exact solution (figure 2) below. Since the exact solution has exponential decay I would expect to see that in the approximation, but it just oscillates between -1 and 1 like a regular sine function. Can anyone point me towards what is wrong with my program? I think it has to do with my choices for n, m, and tf but I’m not sure. How can I fix this without changing the solution method? Any advice is appreciated!
n = 200; m = 100; tf = .5;
x = linspace(0,1,n+1)’; xs = x(2:end-1); dx = 1/n;
t = linspace(0,tf,m+1)’; dt = tf/m;
[X,Y] = meshgrid(x,x); [Xs,Ys] = meshgrid(xs,xs);
% Building the Laplacian
d = kron(ones(n,1),[1,-2,1]);
D2 = (1/dx^2) * spdiags(d, -1:1, n-1, n-1);
I = speye(n-1);
L = kron(I,D2) + kron(D2,I);
% Setting up function
F = @(t,w) [w((n-1)^2+1:end); L*w(1:(n-1)^2)];
% Initial conditions
f = @(x,y) sin(pi*x) .* sin(pi*y);
w0 = [reshape(f(Xs,Ys),[(n-1)^2,1]); zeros((n-1)^2,1)];
% Runge-Kutta
w = RK_solver(F,t,w0);
% Plotting
figure(1);
for i = 1:m+1
colormap(spring);
u = reshape(w(1:(n-1)^2,i),n-1,n-1);
z = zeros(n+1);
z(2:end-1,2:end-1) = u;
surf(X,Y,z); zlim([-1 1]);
pause(0.05)
end
% Exact for comparison
figure(2);
for i = 0:dt:tf
colormap(autumn);
U = exp(-2*pi^2*i) .* sin(pi*X) .* sin(pi*Y);
surf(X,Y,U); zlim([0 1]);
pause(0.05);
end
And here’s the Runge-Kutta solver I used:
function w = RK_solver(F,t,c)
n = length(t)-1; % number of subdivisions
m = length(c); % number of equations
dt = t(2)-t(1); % change in t
w = zeros(m,n+1); % to store answers
w(:,1) = c; % initial condition
% RK loop
for i = 1:n
k1 = dt.*F(t(i), w(:,i));
k2 = dt.*F(t(i)+(dt/2), w(:,i)+(k1/2));
k3 = dt.*F(t(i)+(dt/2), w(:,i)+(k2/2));
k4 = dt.*F(t(i)+dt, w(:,i)+k3);
% Weighted average of k values
k = (k1+(2*k2)+(2*k3)+k4)/6;
% Fills in next column of w
w(:,i+1) = w(:,i) + k;
endHello, I am trying to create a wave equation solver using the method of lines and Runge-Kutta to numerically approximate the solution. If my initial conditions are u(0,x,y)=sin(pix)sin(piy), u_t(0,x,y)=0, then I know the exact solution to be e^(-2pi^2)sin(pix)sin(piy). I plot both my approximation (figure 1) and the exact solution (figure 2) below. Since the exact solution has exponential decay I would expect to see that in the approximation, but it just oscillates between -1 and 1 like a regular sine function. Can anyone point me towards what is wrong with my program? I think it has to do with my choices for n, m, and tf but I’m not sure. How can I fix this without changing the solution method? Any advice is appreciated!
n = 200; m = 100; tf = .5;
x = linspace(0,1,n+1)’; xs = x(2:end-1); dx = 1/n;
t = linspace(0,tf,m+1)’; dt = tf/m;
[X,Y] = meshgrid(x,x); [Xs,Ys] = meshgrid(xs,xs);
% Building the Laplacian
d = kron(ones(n,1),[1,-2,1]);
D2 = (1/dx^2) * spdiags(d, -1:1, n-1, n-1);
I = speye(n-1);
L = kron(I,D2) + kron(D2,I);
% Setting up function
F = @(t,w) [w((n-1)^2+1:end); L*w(1:(n-1)^2)];
% Initial conditions
f = @(x,y) sin(pi*x) .* sin(pi*y);
w0 = [reshape(f(Xs,Ys),[(n-1)^2,1]); zeros((n-1)^2,1)];
% Runge-Kutta
w = RK_solver(F,t,w0);
% Plotting
figure(1);
for i = 1:m+1
colormap(spring);
u = reshape(w(1:(n-1)^2,i),n-1,n-1);
z = zeros(n+1);
z(2:end-1,2:end-1) = u;
surf(X,Y,z); zlim([-1 1]);
pause(0.05)
end
% Exact for comparison
figure(2);
for i = 0:dt:tf
colormap(autumn);
U = exp(-2*pi^2*i) .* sin(pi*X) .* sin(pi*Y);
surf(X,Y,U); zlim([0 1]);
pause(0.05);
end
And here’s the Runge-Kutta solver I used:
function w = RK_solver(F,t,c)
n = length(t)-1; % number of subdivisions
m = length(c); % number of equations
dt = t(2)-t(1); % change in t
w = zeros(m,n+1); % to store answers
w(:,1) = c; % initial condition
% RK loop
for i = 1:n
k1 = dt.*F(t(i), w(:,i));
k2 = dt.*F(t(i)+(dt/2), w(:,i)+(k1/2));
k3 = dt.*F(t(i)+(dt/2), w(:,i)+(k2/2));
k4 = dt.*F(t(i)+dt, w(:,i)+k3);
% Weighted average of k values
k = (k1+(2*k2)+(2*k3)+k4)/6;
% Fills in next column of w
w(:,i+1) = w(:,i) + k;
end Hello, I am trying to create a wave equation solver using the method of lines and Runge-Kutta to numerically approximate the solution. If my initial conditions are u(0,x,y)=sin(pix)sin(piy), u_t(0,x,y)=0, then I know the exact solution to be e^(-2pi^2)sin(pix)sin(piy). I plot both my approximation (figure 1) and the exact solution (figure 2) below. Since the exact solution has exponential decay I would expect to see that in the approximation, but it just oscillates between -1 and 1 like a regular sine function. Can anyone point me towards what is wrong with my program? I think it has to do with my choices for n, m, and tf but I’m not sure. How can I fix this without changing the solution method? Any advice is appreciated!
n = 200; m = 100; tf = .5;
x = linspace(0,1,n+1)’; xs = x(2:end-1); dx = 1/n;
t = linspace(0,tf,m+1)’; dt = tf/m;
[X,Y] = meshgrid(x,x); [Xs,Ys] = meshgrid(xs,xs);
% Building the Laplacian
d = kron(ones(n,1),[1,-2,1]);
D2 = (1/dx^2) * spdiags(d, -1:1, n-1, n-1);
I = speye(n-1);
L = kron(I,D2) + kron(D2,I);
% Setting up function
F = @(t,w) [w((n-1)^2+1:end); L*w(1:(n-1)^2)];
% Initial conditions
f = @(x,y) sin(pi*x) .* sin(pi*y);
w0 = [reshape(f(Xs,Ys),[(n-1)^2,1]); zeros((n-1)^2,1)];
% Runge-Kutta
w = RK_solver(F,t,w0);
% Plotting
figure(1);
for i = 1:m+1
colormap(spring);
u = reshape(w(1:(n-1)^2,i),n-1,n-1);
z = zeros(n+1);
z(2:end-1,2:end-1) = u;
surf(X,Y,z); zlim([-1 1]);
pause(0.05)
end
% Exact for comparison
figure(2);
for i = 0:dt:tf
colormap(autumn);
U = exp(-2*pi^2*i) .* sin(pi*X) .* sin(pi*Y);
surf(X,Y,U); zlim([0 1]);
pause(0.05);
end
And here’s the Runge-Kutta solver I used:
function w = RK_solver(F,t,c)
n = length(t)-1; % number of subdivisions
m = length(c); % number of equations
dt = t(2)-t(1); % change in t
w = zeros(m,n+1); % to store answers
w(:,1) = c; % initial condition
% RK loop
for i = 1:n
k1 = dt.*F(t(i), w(:,i));
k2 = dt.*F(t(i)+(dt/2), w(:,i)+(k1/2));
k3 = dt.*F(t(i)+(dt/2), w(:,i)+(k2/2));
k4 = dt.*F(t(i)+dt, w(:,i)+k3);
% Weighted average of k values
k = (k1+(2*k2)+(2*k3)+k4)/6;
% Fills in next column of w
w(:,i+1) = w(:,i) + k;
end pde, differential equations, numerical analysis, numerical approximation, runge-kutta, runge kutta, method of lines, partial differential equations, wave equation, surf MATLAB Answers — New Questions