Solving the Wave Equation (DD case) using the Discrete Sine Transform

Contents

You first have to download m-files

Download the m-files sintr.m , isintr.m , colorcurves.m and put them in the same directory as your other m-files.

Solve wave equation $u_{tt} = u_{xx}$

We want to find a function $u(x,t)$ for $x\in [0,1]$, $t\in[0,\infty)$ such that

$$u_{tt}(x,t) = u_{xx}(x,t) \qquad x\in (0,1)$, $t\in(0,\infty) \quad\mbox{(PDE)}$$

$$u(x,0) = f(x), \quad u_t(x,0) = g(x) \qquad x\in (0,1) \quad\mbox{(Initial Conditions)}$$

where $f,g$ are given function on $(0,1)$.

Boundary conditions: At the endpoints we have DIRICHLET conditions $u(0,t)=0$ and $u(1,t)=0$.

The method of Separation of Variables gives the following Solution Formula:

First find the sine series $f(x)=\sum_{k=1}^\infty b_k \sin(k\pi x)$ and $g(x)=\sum_{k=1}^\infty c_k \sin(k\pi x)$. Let $\mu_k:=k\pi$ for $k=1,2,\ldots$. Then the solution $u(x,t)$ is given by

$$u(x,t) = \sum_{k=1}^\infty \sin(\mu_k x) \left[ b_k \cos(\mu_k t) + \frac{c_k}{\mu_k}\sin(\mu_k t) \right]$$

Example: given initial functions $f(x)= \max\{.1-|x-.3|,0\}$ and $g(x)=0$.

We see that the initial displacement causes waves going to the left and right. These waves are reflected at the boundaries. At time $t=1$ the solution looks similar to $t=0$, but instead of a positive peak at .3 we have a negative peak at .7.

f = @(x) max(.1-abs(x-.3),0);  % given initial function f
N = 100;                       % use spacing 1/100 for x
x = (1:N-1)'/N;                % only interior nodes            (column vector)
F = f(x);                      % evaluate f for 1/N,...,(N-1)/N (column vector)
plot(x,F); title('initial condition u(x,0)')

Find solution using Discrete Sine Transform

b = sintr(F);                  % take Discrete Sine Transform

t = 0:.02:1;                   % find solution for these t-values (row vector)
kv = (1:N-1)';                 % column vector of k-values
muv = kv*pi;                   % column vector of mu-values
W = cos(muv*t);                % array of cos(mu_k*t) for k=1...N-1, for all t-values
                               % (column vector)*(row vector) gives array
U = isintr(diag(b)*W);         % multiply 1st row by b(1), 2nd row by b(2), ...
                               % then take Inverse Sine Transform of each column

colorcurves(x,t,U); axis tight
xlabel('x'); ylabel('t'); title('u(x,t)')
view(150,30)

Solution at time $t=1$

plot(x,U(:,end),'r'); title('solution at time t=1')


Published with MATLAB® R2018a