Numerical solution to the wave equation-Effect of grid size

OBJECTIVE- To write a code to solve the 1D linear wave equation equation numerically.

1. Assume that the domain length is L = 1m

2. The initial velocity profile is a step function. It is equal to 2m/s between x= 0.1 and 0.3 and 1m/s everywhere else

3. Use first order forward differencing for the time derivative

4. Use first order backward/rearward differencing for the space term.

5. Time step and number of grid points should be variables in your calculation

6. C = 1

7. Set time step = 0.01

Expected output:

1. For C = 1 and for time = 0.4seconds compare the original and final velocity profiles.

2. Make the comparison for n=20,40,80 and 160

1D linear wave equation

`(delu)/(delt)+c*(delu)/(delx)=0`

Discretizing the equation by using forward differencing scheme(FDS) for time derivative and backward differencing scheme(BDS) for space term we get,

`(u_i^(n+1)-u_i^n)/(Deltat)=-C*(u_i^n-u_(i-1)^n)/(Deltax)`

Plot for velocity vs domain length for n=20

From above graph we can see that the velocity profile is moving along the length of domain.The peak of the profile is dying.The reason for peak getting reduced is the error in the numerical scheme which is causing artificial diffusion.The final velocity profile is drawn at 0.4 sec.

Plot for velocity vs domain length for n=40

As we have increased the no of grids peak of the profile is very slowly dying.The final velocity profile is drawn at 0.4 sec.

Plot for velocity vs domain length for n=80

For the above graph the no of grid points is 80.And we can see that the peak of profile is maintained as it moves along the domain.And the intial velocity profile is getting more squarer.

Plot for velocity vs domain length for n=160

when n becomes 160 the solution seems to be blown off.It indicates that the numerical sheme is not stable.

clear all;
close all;
clc;

L=1;
n=160;
c=1;

nstart=0;nstop=0;
x=linspace(0,1,n);
      
 dx=x(2)-x(1);
 nstart=(0.1/dx)+1;
  nstop=(0.3/dx)+1; 
dt=0.01;
t_end=40;
u=ones(1,n);

u(nstart:nstop)=2;
u(1)=1;
 uold=u;

 uinitial=u;

for k=1:t_end;
   
for j=2:length(x);
    u(j)=uold(j)-c*(dt/dx)*(uold(j)-uold(j-1));
end

uold=u;
grid on
plot(x,u,'r');
hold on
plot(x,uinitial,'g','linewidth',2);
 axis([0 1.2 0 2.5]);
 pause(0.03);
 xlabel('Domain length')
 ylabel('velocity')
 title('n=160')
 legend('Final velocity after 0.4sec','initial velocity')
end