- Use mesh/meshgrid to define the u function in x and t.
- Re-map t: t=linspace(0,1,101), then t(30)=0.3
Solving a wave equation in matlab
56 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Here is the problem statement:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/283202/image.png)
I am having trouble plotting the solution at t=0.3 and not sure if my code is solving the wave equation correctly. This is what I have, but not sure where to go from here. Any help would be great thanks!
xstep = 0.1;
tstep = 0.05;
xstep2 = xstep*xstep;
tstep2 = tstep*tstep;
alpha = 2;
alpha2 = alpha*alpha;
lambda2 = alpha2*tstep2/xstep2;
xdomain = [0 1];
tdomain = [0 1];
nx = round((xdomain(2)-xdomain(1))/xstep);
nt = round((tdomain(2)-tdomain(1))/tstep);
xt0 = zeros((nx+1),1); % initial condition
dxdt0 = zeros((nx+1),1); % initial derivative
xold = zeros((nx+1),1); % solution at timestep k
x2old = zeros((nx+1),1); % solution at timestep k-1
xnew = zeros((nx+1),1); % solution at timestep k+1
% initial condition
pi = acos(-1.0);
for i=1:nx+1
xi = (i-1)*xstep;
if(xi>=0 && xi<=1)
xt0(i) = sin(2*pi*xi);
dxdt0(i) = alpha*pi*sin(2*pi*xi);
xold(i) = xt0(i)+dxdt0(i)*tstep;
xold(i) = xold(i) - 4*pi*pi*sin(2*pi*xi)*tstep2*alpha2;
end
end
close all
syms x
t=0.3;
x=linspace(xdomain(1),xdomain(2),nx+1);
analy= sin(2*pi*x)*(sin(4*pi*t)+cos(4*pi*t));
h1=plot(x,analy,'linewidth',2);
hold on;
h2=plot(x,xold(:,t),'linewidth',2);
hold on;
h3=plot(x,xnew(:,t),'linewidth',2);
hold off
legend('Analytical','Initial','Final')
xlabel('x [m]');
ylabel('Displacement [m]');
set(gca,'FontSize',16);
for k=2:nt
time = i*tstep;
for i=1:nx+1
% Use periodic boundary condition, u(nx+1)=u(1)
if(i==1)
xnew(i) = 2*(1-lambda2)*xold(i) + lambda2*(xold(i+1)+xold(nx+1)) - x2old(i);
elseif(i==nx+1)
xnew(i) = 2*(1-lambda2)*xold(i) + lambda2*(xold(1)+xold(i-1)) - x2old(i);
else
xnew(i) = 2*(1-lambda2)*xold(i) + lambda2*(xold(i+1)+xold(i-1)) - x2old(i);
end
end
x2old=xold;
xold = xnew;
if(mod(k,2)==0)
h3.YData = xnew;
refreshdata(h3);
pause(0.5);
end
end
1 comentario
Optics Wizard
el 9 de Abr. de 2020
I haven't gone through your code on the full-scale, but your time-mapping is currently incorrect. You can only index values by integers, but you're currently indexing to data point 0.3:
h2=plot(x,xold(:,t),'linewidth',2);
For instance, when t=0.3, this line returns an error.
Some options:
I hope that helps!
Respuestas (1)
Sulaymon Eshkabilov
el 30 de Oct. de 2021
Here is a part of your code that is corrected:
xstep = 0.1;
tstep = 0.05;
xstep2 = xstep*xstep;
tstep2 = tstep*tstep;
alpha = 2;
alpha2 = alpha*alpha;
lambda2 = alpha2*tstep2/xstep2;
xdomain = [0 1];
tdomain = [0 1];
nx = round((xdomain(2)-xdomain(1))/xstep);
nt = round((tdomain(2)-tdomain(1))/tstep);
xt0 = zeros((nx+1),1); % initial condition
dxdt0 = zeros((nx+1),1); % initial derivative
xold = zeros((nx+1),1); % solution at timestep k
x2old = zeros((nx+1),1); % solution at timestep k-1
xnew = zeros((nx+1),1); % solution at timestep k+1
% initial condition
pi = acos(-1.0);
for i=1:nx+1
xi = (i-1)*xstep;
if(xi>=0 && xi<=1)
xt0(i) = sin(2*pi*xi);
dxdt0(i) = alpha*pi*sin(2*pi*xi);
xold(i) = xt0(i)+dxdt0(i)*tstep;
xold(i) = xold(i) - 4*pi*pi*sin(2*pi*xi)*tstep2*alpha2;
end
end
close all
syms x
t=0.3;
x=linspace(xdomain(1),xdomain(2),nx+1);
analy= sin(2*pi*x)*(sin(4*pi*t)+cos(4*pi*t));
h1=plot(x,analy,'linewidth',2);
hold on;
h2=plot(x,xold(:,end),'linewidth',2); % Index issue in xold
hold on;
h3=plot(x,xnew(:,end),'linewidth',2); % Index issue in xnew
hold off
legend('Analytical','Initial','Final')
xlabel('x [m]');
ylabel('Displacement [m]');
set(gca,'FontSize',16);
for k=2:nt
time = i*tstep;
for i=1:nx+1
% Use periodic boundary condition, u(nx+1)=u(1)
if(i==1)
xnew(i) = 2*(1-lambda2)*xold(i) + lambda2*(xold(i+1)+xold(nx+1)) - x2old(i);
elseif(i==nx+1)
xnew(i) = 2*(1-lambda2)*xold(i) + lambda2*(xold(1)+xold(i-1)) - x2old(i);
else
xnew(i) = 2*(1-lambda2)*xold(i) + lambda2*(xold(i+1)+xold(i-1)) - x2old(i);
end
end
x2old=xold;
xold = xnew;
if(mod(k,2)==0)
h3.YData = xnew;
refreshdata(h3);
pause(0.5);
end
end
2 comentarios
Ver también
Categorías
Más información sobre Eigenvalue Problems en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!