Plotting wave solution at specific time

Question

bml727 2020-4-16

0
链接

此问题的直接链接

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

回答： Deepak 2024-11-13，6:34

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Deepak 2024-11-13，6:34

0
链接

此回答的直接链接

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

在 MATLAB Online 中打开

Hi @bml727,

We can plot the solution at time step (t = 0.3) by updating the plot for the sixth iteration of the solution array.

We can update the plot when the loop reaches this specific time step, so that the “xnew” values are correctly displayed. This involves setting the “YData” of the plot to the solution at this iteration, ensuring “xnew” values are correctly displayed on graph and adding a title to indicate the time step.

Here is the updated MATLAB code with the required changes:

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); 
analy = sin(2*pi.*x).*(sin(4*pi.*0.3) + cos(4*pi.*0.3)); 
h1 = plot(x, analy, 'DisplayName', 'Analytical'); 
hold on; 
h2 = plot(x, xt0, 'linewidth', 2, 'DisplayName', 'Initial'); 
hold on; 
h3 = plot(x, xnew, 'DisplayName', 'Final'); 
legend('show'); 
xlabel('x [m]'); 
ylabel('Displacement [m]'); 
set(gca, 'FontSize', 16); 
tplot = zeros(1, nt); 
for k = 1:nt 
    time = k * tstep; 
    tplot(k) = time; 
    
    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; 
    
    % Update plot at t = 0.3 (tplot(6)) 
    if k == 6 
        h3.YData = xnew; 
        title('Solution at t = 0.3'); 
        break; % Exit the loop after updating the plot for t = 0.3 
    end 
end 