How to plot a simple wave equation using method of lines?

Question

Iris 2024-4-14，6:08

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2106586-how-to-plot-a-simple-wave-equation-using-method-of-lines

编辑： Torsten 2024-4-14，19:44

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

2 个评论
显示无隐藏无

Torsten 2024-4-14，18:07

编辑：Torsten 2024-4-14，18:07

Where do you plot the red solution curves in your code ? I can only see the exact solution in blue over time.

Iris 2024-4-14，18:24

在 MATLAB Online 中打开

Sorry, I left that out on accident. It was just this loop:

for i = 1:m+1
    u = w(:,i);
    plot(x,[0;u;0],'r'); hold on 
    pause(0.01);
end

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-4-14，19:04

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2106586-how-to-plot-a-simple-wave-equation-using-method-of-lines#answer_1441016

编辑：Torsten 2024-4-14，19:31

在 MATLAB Online 中打开

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

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

Torsten 2024-4-14，19:24

编辑：Torsten 2024-4-14，19:44

My guess is there is something wrong with the system of differential equations you pass to the Runge-Kutta function. The integrator looks correct (see above for a solution with your Runge-Kutta code, but my ODE system).

Iris 2024-4-14，19:29