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

Answer 1

[Deepak]

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 

Attached is the documentation of functions referenced:

plot: www.mathworks.com/help/matlab/ref/plot.html

acos: www.mathworks.com/help/matlab/ref/double.acos.html

linspace: www.mathworks.com/help/matlab/ref/double.linspace.html

I hope this helps in resolving the issue.