Solving a wave equation in matlab

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

[Optics Wizard]

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:

1. Use mesh/meshgrid to define the u function in x and t.
2. Re-map t: t=linspace(0,1,101), then t(30)=0.3

I hope that helps!

Answer 1

[Sulaymon Eshkabilov]

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

[Samia Alaa]

Can you helpe me to solve part one find A and B

[naren BORO]

final part of the graph is not coming. please help me.