$$
{\bf u}^{n+1} = {\bf u}^{n} – \frac{\Delta t}{2 \Delta x} {\bf c}.^*({\bf D}_{\bf x}{\bf u}^n) + \frac{\Delta t^2}{2 \Delta x^2} {\bf c}^2.^*({\bf D}_{\bf x x}{\bf u}^n) + \frac{\Delta t^2}{8 \Delta x^2} {\bf c}.^*({\bf D}_{\bf x}{\bf c}).^*({\bf D}_{\bf x}{\bf u}^n) ,
$$
where ${\bf c}^2$ is the element-wise square of the vector ${\bf c}$. You should state clearly what the vectors ${\bf u}^n$ and ${\bf c}$ and matrices ${\bf D}_{\bf x}$ and ${\bf D}_{\bf xx}$ are.Take $x \in [0,10]$, $u(0,t) = u(10,t) = 0$, $c = 1/5$ and initial condition $u(x,0) = e^{-100 (x-5)^2}$. Using MATLAB, numerically implement the Lax-Wendroff scheme for $N = 50$ and plot the solution at $t = 1$. You are required to choose a time step such that the scheme is numerically stable. Discuss your reasoning for the time step you chose.
I am trying to code a solution for this PDE $ u_t +c(x)u_x = 0 $ in Matlab using the above scheme, However I am really unsure about how to define my boundary conditions in Matlab. Since the boundary equals zero, but not sure if i start at $j=1$ or $j = 0$, which would mean i have to introduce ghost pints, which i wasn't been given, so i am thinking i don't need to this?
Also Matlab says my matrix dimensions dont agree, but i cant seem to fix this? Can anyone see the error?
My code so far,
Nx = 100;
x = linspace(0,10,Nx);
dx = 10/Nx;
tmax = 1;
dt = dx/50; % ensure stability
R = round(tmax/dt);
e = ones(Nx,1);
Dx = spdiags([-e e],[-1 1],Nx,Nx);
Dxx = spdiags([e,-2*e,e],[-1 1], Nx, Nx);
Dx(1,2) = 0;
Dx(Nx,Nx-1) = 0; % change entries determined by BCs
c = 1/5*ones(Nx,1);
p = 0.2*dt/dx;
% initial data
u = zeros(Nx,1);
% initial condition
for i=1:11
u(:,1,i)=exp(-100*(x-0.5).^2)';
end
for n = 2:R
u(:,n+1)=u(:,n) - 1\2*p*c.*Dx.*u(:,u) + 1/2*p^2*c^2.*Dxx*u(:,u) + 1/8*p^2*c^2.*Dx*u(:,u);
end
plot(x,u(:,),'-or');
xlabel('x');
ylabel('u');
title('Solution at t = 1');
Best Answer
We refer to this post for notations. The spatial discretization is $x_j = j \Delta x$ for $j = 0,\dots , N$, where $\Delta x = 10/N$. Since $u_0^n \approx u(0,t)$ and $u_N^n \approx u(10,t)$ are known, the vector of unknowns ${\bf u}^n = (u_1^n , \dots , u_{N-1}^n)^\top$ is of size $N-1$. Moreover, since $u_0^n$ and $u_N^n$ equal zero, they do not appear in the right-hand side of the scheme at the lines corresponding to the computation of $u_1^{n+1}$ and $u_{N-1}^{n+1}$. Lastly, let us introduce the vector ${\bf c} = (c_1, \dots , c_{N-1})^\top$. Thus, using MATLAB's matrix notation, $$ {\bf u}^{n+1} = {\bf u}^{n} - \frac{\Delta t}{2 \Delta x} {\bf c}.^*({\bf D}_{\bf x}{\bf u}^n) + \frac{\Delta t^2}{2 \Delta x^2} {\bf c}^2.^*({\bf D}_{\bf x x}{\bf u}^n) + \frac{\Delta t^2}{8 \Delta x^2} {\bf c}.^*({\bf D}_{\bf x}{\bf c}).^*({\bf D}_{\bf x}{\bf u}^n) , $$ where ${\bf D}_{\bf x}$ and ${\bf D}_{\bf xx}$ are the $({N-1})\times ({N-1})$ matrices given by $$ {\bf D}_{\bf x} = \left( \begin{array}{cccc} 0 & 1 & & \\ -1 & \ddots & \ddots & \\ & \ddots & \ddots & 1 \\ & & -1 & 0 \\ \end{array} \right) , \quad {\bf D}_{\bf xx} = \left( \begin{array}{cccc} -2 & 1 & & \\ 1 & \ddots & \ddots & \\ & \ddots & \ddots & 1 \\ & & 1 & -2 \\ \end{array} \right) . $$ Here is a corresponding code in MATLAB syntax:
Since the speed $c = 1/5$ is constant, the result at $t=1$ should be a translation of the initial data by $ct \approx 0.2$, i.e. $u(x,t) = u(x-ct,0)$. Despite the fact that $N=50$ is not large enough for a good resolution, this analytical result (black solid line) is coherent with the numerical result (blue circles):
Here, the Courant number is set to $0.95$. As shows von Neumann analysis, stability is obtained for Courant numbers smaller than one. The Courant number should be as large as possible so as to increase $\Delta t$ and reduce the number of iterations in time. Moreover, numerical diffusion increases when the Courant number diminishes.