How to plot a simple wave equation using method of lines?
56 views (last 30 days)
Show older comments
Hello, I am trying to build a simple wave equation solver that uses the method of lines and runge-kutta. The example I am testing is of the form u_tt=u_xx with initial conditions u(0,x)=sin(2pix), u_t(0,x) = 0. I know the exact solution to be cos(2pi*t)sin(2pi*x), whose range is from -1 to 1. However, running my solver gives me very large values, and I'm not sure where they are coming from. Here is my approximation (red), vs the exact solution (blue):
I think these values are way too big, even considering error. Can anyone point me to where I've gone wrong? I am not super familiar with the wave equation. Here is my program and Runge-Kutta solver:
n = 100; m = 100; tf = 1;
x = linspace(0,1,n+1)'; xs = x(2:end-1); dx = 1/n; dt = 1/m;
f = @(x) sin((2*pi).*x); g = 0; % initial conditions
exact = @(t,x) (cos((2*pi).*t))*(sin((2*pi).*x)); % exact for comparison
% building second difference matrix
d = kron(ones(n,1),[1,-2,1]); D2 = (1/dt^2) * spdiags(d, -1:1, n-1, n-1);
% build system of eqns
w0 = [f(xs);zeros(size(xs))];
t = linspace(0,tf,m+1)'; F = @(x,w) [w(n:end);D2*w(1:n-1)];
w = RK_solver(F,t,w0); w = reshape(w,n-1,(m+1)*2);
figure(1)
hold on
for i = 1:m+1
%u = w(:,i);
plot(x,exact(t(i),x),'b'); axis([0 1 -1 1])
%pause(0.01);
end
hold off
figure(2)
hold on
for i = 1:m+1
u = w(:,i);
plot(x,[0;u;0],'r'); %hold on
%pause(0.01);
end
hold off
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);
w(:,1) = c; % To store answers
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);
k = (k1+(2*k2)+(2*k3)+k4)/6;
w(:,i+1) = w(:,i) + k;
end
end
2 Comments
Accepted Answer
Torsten
on 14 Apr 2024 at 19:04
Edited: Torsten
on 14 Apr 2024 at 19:31
n = 100; m = 100; tf = 1; dx = 1/n;
x = linspace(0,1,n+1)';
tspan = linspace(0,tf,m+1);
y0 = [sin(2*pi*x);zeros(n+1,1)];
%[t,y] = ode15s(@(t,y)fun(t,y,n,dx),tspan,y0);
y = RK_solver(@(t,y)fun(t,y,n,dx),tspan,y0);
y = y';
exact = cos(2*pi*tspan.').*sin(2*pi*x.');
figure(1)
hold on
for i = 1:m+1
plot(x,exact(i,:),'b')
end
hold off
axis([0 1 -1 1])
figure(2)
hold on
for i = 1:m+1
plot(x,y(i,1:n+1),'r')
end
hold off
axis([0 1 -1 1])
function dy = fun(t,y,n,dx)
u = y(1:n+1);
v = y(n+2:2*(n+1));
dy = zeros(2*(n+1),1);
dy(1:n+1) = v;
dy(n+2:2*(n+1)) = [0;(u(3:n+1)-2*u(2:n)+u(1:n-1))/dx^2;0];
end
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);
w(:,1) = c; % To store answers
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);
k = (k1+(2*k2)+(2*k3)+k4)/6;
w(:,i+1) = w(:,i) + k;
end
end
3 Comments
More Answers (0)
See Also
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!