Plotting wave solution at specific time

1 次查看（过去 30 天）

bml727 2020-4-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/518434-plotting-wave-solution-at-specific-time

I solved an equation u(x,t) and need to plot the solution at t=0.3, but I am not sure how to do it. I saved the time iteration and found that t=0.3 occurs at tplot(6). Any help would be great!

clear;
clc;
%% Problem 2 
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 
x2old = xt0;
close all
x=linspace(xdomain(1),xdomain(2),nx+1);
%t=linspace(tdomain(1),tdomain(2),nx+1);
analy= sin(2*pi.*x).*(sin(4*pi.*0.3+cos(4*pi.*0.3)));
h1=plot(x,analy);
hold on;
h2=plot(x,xt0,'linewidth',2);
hold on;
h3=plot(x,xnew);
legend('Analytical','Initial','Final')
xlabel('x [m]');
ylabel('Displacement [m]');
set(gca,'FontSize',16);
tplot=zeros(1,nt);
for k=1:nt
    time = k*tstep;
     tplot(k)=time;
     %tplot(6)=0.3
    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