- 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
66 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Here is the problem statement:
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 commentaire
Optics Wizard
le 9 Avr 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!
Réponses (1)
Sulaymon Eshkabilov
le 30 Oct 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 commentaires
Voir également
Catégories
En savoir plus sur Eigenvalue Problems dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!