Solving a wave equation in matlab

Question

bml727 2020-4-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/516731-solving-a-wave-equation-in-matlab

评论： naren BORO 2022-6-2

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Optics Wizard 2020-4-9

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:

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

I hope that helps!

请先登录，再进行评论。

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

Answer 1

Sulaymon Eshkabilov 2021-10-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/516731-solving-a-wave-equation-in-matlab#answer_820143

在 MATLAB Online 中打开

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

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

x2old=xold;

xold = xnew;

if(mod(k,2)==0)

h3.YData = xnew;

refreshdata(h3);

pause(0.5);

end

2 个评论
显示无隐藏无

Samia Alaa 2021-10-30

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

naren BORO 2022-6-2

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

请先登录，再进行评论。

Solving a wave equation in matlab

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

回答（1 个）

2 个评论
显示无隐藏无

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Solving a wave equation in matlab

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

回答（1 个）

2 个评论 显示 无隐藏 无

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

WeChat

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

2 个评论
显示无隐藏无