1D wave equation using Finite difference method MATLAB

I have the wave equation $$u_{tt} = 4 u_{xx}$$ with the boundary conditions $$u(0,t) = u(L,t) = 0\,,\quad x \leq 0 \leq 2\pi \,,\quad t\geq 0$$ and initial conditions \begin{align} &u(x,0)=\begin{cases}1 & \text{for} \quad \pi-1 \leq x \leq \pi +1 \\0 &\text{otherwise} \end{cases} \\[0.5em] &u_{t}(x,0) = x^2sin^2(4x) \end{align}

I'm trying to implement this problem on MATLAB by the finite difference method and by using the surf function to plot it as a 3D wave; however, the problem I'm having is how to code the first initial condition.

The function I am using is

function [u, q] = Wave(f1,f2,g0,g1,xspan,tspan,nx,nt,alpha)
x0 = xspan(1)
xf = xspan(2)
t0 = tspan(1)
tf = tspan(2)
dx = (xf - x0)/nx;
dt = (tf-t0)/nt;
x = [0:nx]'*dx;
t = [0:nt]*dt;
q = alpha*(dt/dx)^2;
q1 = q/2;
q2 = 2*(1-q);
u(:,1) = f1(x);
for k = 1:nt+1, u([1 nx+1],k) = [g0(t(k)); g1(t(k))]; end
u(2:nx,2) = q1*u(1:nx-1,1) + (1-q)*u(2:nx,1) + q1*u(3:nx+1,1) + dt*f2(x(2:nx));
for k = 3:nt+1
u(2:nx,k) = q*u(1:nx-1,k-1) + q2*u(2:nx,k-1) + q*u(3:nx+1,k-1) - u(2:nx,k-2);
end

surf(t,x,u)
xlabel('t')
ylabel('x')
zlabel('u(x,t)')
end

and I know that within the function, the variable f1 is the one controlled by the first initial condition. In the function, I don't know how to incorporate the initial condition into the method. I assume that a for loop and an if statement are to be used but anything I tried doesn't work.

The script that I am using to plot the graph is

clc
clear all
f1 = @(x) ////;
f2 = @(x) (x.^2).*((sin(4.*x)).^2);
g0 = @(t) 0;
g1 = @(t) 0;
xspan = [0,2*pi];
tspan = [0,1];
nx = 20;
nt = 40;
alpha = 4;
[u,r] = Wave(f1,f2,g0,g1,xspan,tspan,nx,nt,alpha);

As seen on the script where f1 = @(x) ////;. I'm unsure of what to put here for the initial condition. Any help would be very much appreciated.

• Well you need to set u(:, 1) to something otherwise your time stepping won’t have something to initially refer to. Nov 20, 2021 at 22:45

Check the interval of x and then covert logical to double.

f1 = @(x) double(pi-1 <= x && x <= pi+1);

Or if the initial condition is not 1 and 0, write a function

function u = f1(x)
u = zeros(size(x));
is_one = pi-1 <= x && x <= pi+1;
u(is_one) = 1;
% u(~is_one) = 0;  not required, just for your reference
end

Edit:

$$u(x, 0)$$ is set by

u(1,:) = f1(x);

Then you used it to calculate $$u(x,\Delta t)$$, i.e. u(:,2) in this line

u(2:nx,2) = q1*u(1:nx-1,1) + (1-q)*u(2:nx,1) + q1*u(3:nx+1,1)...
+ dt*f2(x(2:nx));

You have done it, haven't you?

• I added u(1,:) = f1(x) to the function and I added f1 = @(x) double(pi-1 <= x && x <= pi+1); to the script but I'm not sure how to implement the function to the program. Nov 21, 2021 at 14:17
• do I change the f2 to f1 on the line highlighted in the edit such as u(2:nx,2) = q1*u(1:nx-1,1) + (1-q)*u(2:nx,1) + q1*u(3:nx+1,1)... + dt*f1(x(2:nx)); ? Nov 22, 2021 at 18:28
• No. Compare this to the formulation of finite difference method. dt*f1(x(2:nx)) means $\Delta t \cdot u_t(x, 0)$, this is integration with "finite difference". If you change f2 to f1 it becomes $\Delta t \cdot u(x, 0)$, doesn't make much sense. Nov 23, 2021 at 1:09
• Do I have to add f1 into it somehow because the plotted graph doesn’t look correct even with the changes. Nov 24, 2021 at 3:42
• You have to turn back to the algorithm of finite difference method and check if your code match the equations in the algorithm. Pay attention to how it evloves from this time step to the next one. Nov 25, 2021 at 9:20