Aproximación de la ecuación de ondas con condiciones Dirichlet: Método de líneas

This article shows how to approximate the wave equation with the method of lines.

1 Wave equation

Wave equation is used to simulate a number of applications related with waves, as the vibration of strings. In the one dimensional case it reads,

[math]u_{tt}(x,t) -u_{xx}(x,t)= f(x,t), \quad x\in(a,b), \quad t\in(0,T) [/math]

[math]u(a,t) = g_1(t), \qquad u(b,t)=g_2(t), \quad t\in(0,T) [/math]

[math]u(x,0) = u^0(x), \quad x\in(a,b),[/math]

[math]u_t(x,0) = u^1(x), \quad x\in(a,b),[/math]

Here, x is the space variable and t is the time, [math]u(x,t)[/math] represents the vertical displacement of a string located in the segment [math](a,b)[/math], [math]g_1(t),g_2(t)[/math] represent the vertical displacement at the extremes (which is usually a datum) and [math]u^0(x), u^1(x)[/math] the initial configuration and velocity of the points of the string at time [math]t=0[/math].

2 Numerical scheme

Space discretization

Take a partition of [0,L] with step h. Define [math] u_n(t)\sim u(x_n,t),[/math] and replace

[math]-u_{xx}(x,t) \quad \mbox{ by }\quad \frac{-u_{n-1}(t)+2u_n(t) - u_{n+1}(t)}{h^2}[/math]

Then, we have the system

[math]\left\{ \begin{array}{l} u_n''(t)+\frac{-u_{n-1}(t)+2u_n(t) - u_{n+1}(t)}{h^2}=f(x_n,t), \quad n=1,2,...,N-1 \qquad t\gt0, \\ u_0(t)=g_1(t), \qquad t\gt0, \\ u_{N}(t)=g_2(t), \qquad t\gt0, \\ u_n(0)=u^0(x_n), \qquad n=1,2,...,N-1,\\ u_n'(0)=u^1(x_n), \qquad n=1,2,...,N-1. \end{array} \right. [/math]

The first N+1 equations can be written as

[math] \left\{ \begin{array}{l} {u_0(t)=g_1(t)},\\ u_1''(t)+\frac{{-u_0(t)+}2u_1(t) - u_{2}(t)}{h^2}=f(x_1,t), \\ u_2''(t)+\frac{-u_{1}(t)+2u_2(t) - u_{3}(t)}{h^2}=f(x_2,t), \\ u_3''(t)+\frac{-u_{2}(t)+2u_3(t) - u_{4}(t)}{h^2}=f(x_3,t), \\ ... \\ u_{N-1}''(t)+\frac{-u_{N-2}(t)+2u_{N-1}(t) {- u_{N}(t)}}{h^2}=f(x_{N-1},t), \\ {u_N(t)=g_2(t)}, \end{array} \right. [/math]

Eliminate the equations for [math]u_0(t)[/math] and [math]u_{N}(t)[/math],

[math] \left\{ \begin{array}{l} u_1''(t)+\frac{2u_1(t) - u_{2}(t)}{h^2}=f(x_1,t)+g_1(t)/h^2, \\ u_2''(t)+\frac{-u_{1}(t)+2u_2(t) - u_{3}(t)}{h^2}=f(x_2,t), \\ u_3''(t)+\frac{-u_{2}(t)+2u_3(t) - u_{4}(t)}{h^2}=f(x_3,t), \\ ... \\ u_{N-1}''(t)+\frac{-u_{N-2}(t)+2u_{N-1}(t)}{h^2}=f(x_{N-1},t)+g_2(t)/h^2, \end{array} \right. [/math]

The initial condition is tranlated into the 2x(N-1) equations:

[math] u_n(0)=u^0(x_n), \qquad u_n'(0)=u^1(x_n), \qquad n=1,2,...,N-1. [/math]

In matrix form:

[math] U''(t)+KU(t)=F(t), \qquad {U(0)=U^0}, \qquad U'(0)=U^1, [/math]

where [math] K=\frac1{h^2}\left( \begin{array}{ccccccc} 2& -1& 0& 0& ...& 0 &0 \\ -1& 2& -1& 0& ...& 0 &0 \\ ...& ...& ...& ...& ...& ...& ...\\ 0& 0& 0& 0& ...& -1& 2 \end{array} \right), \quad U=\left( \begin{array}{l} u_1\\ u_2 \\... \\u_{N-1} \end{array} \right) [/math]

[math] U^0=\left( \begin{array}{l} u^0(x_1) \\ u^0(x_2)\\ ... \\ u^0(x_{N-1}) \end{array} \right), \quad U^1=\left( \begin{array}{l} u^1(x_1) \\ u^1(x_2)\\ ... \\ u^1(x_{N-1}) \end{array} \right), \qquad F=\left( \begin{array}{l} f(x_1,t)+g_1(t)/h^2 \\ f(x_2,t)\\ ... \\ f(x_{N-1},t)+g_2(t)/h^2 \end{array} \right) [/math]

Time discretization

[math] \left\{ \begin{array}{l} U(0)=U^0, \quad U'(0)=U^1, \\ U''(t)+KU(t)=F(t), \end{array} \right. [/math]

We reduce the system to a first order one by introducing a new variable [math]V=U'[/math]

[math] \left\{ \begin{array}{l} \left( \begin{array}{l} U \\ V \end{array} \right) (0) = \left( \begin{array}{l} U^0 \\ U^1 \end{array} \right), \\ \left( \begin{array}{l} U \\ V \end{array} \right)'(t)+\left( \begin{array}{ll} 0 & -I \\ K & 0 \end{array} \right) \left( \begin{array}{l} U \\ V \end{array} \right) = \left( \begin{array}{l} 0 \\ F(t) \end{array} \right) , \end{array} \right. [/math]

or

[math] \left\{ \begin{array}{l} W (0) = W^0 , \\ W'(t)+AW=G , \end{array} \right. [/math]

where,

[math] W=\left( \begin{array}{l} U \\ V \end{array} \right), \quad W^0=\left( \begin{array}{l} U^0 \\ V^0 \end{array} \right), \quad A=\left( \begin{array}{ll} 0 & -I \\ K & 0 \end{array} \right), \quad G=\left( \begin{array}{l} 0 \\ F(t) \end{array} \right). [/math]

Chose a time mesh [math]\{t^j \}_{j=0}^J[/math] with time step [math]\Delta t[/math] and define

[math] W^j\sim W(t^j), \qquad f_u(t,W)=-A W+G(t). [/math]

Now we replace the ODE by any numerical method for initial value problems. For example, the Euler method: [math] \left\{ \begin{array}{l} W^0, \\ W^{j+1}= W^j+ \Delta t \; f_u(t^j,W^j) , \qquad j=0,1,...,J-1. \end{array} \right. [/math]

Remark Euler method requires [math]\Delta t \leq h[/math] for stability.

Other possible methods are:

1. Implicit Euler (unconditionally stable),

[math] \left\{ \begin{array}{l} W^0, \\ W^{j+1}= W^j+ \Delta t \; f_u(t^{j+1},W^{j+1}), \qquad j=0,1,...,J-1 \end{array} \right. [/math]

2. Trapezoidal (or Crank-Nicolson method, unconditionally stable),

[math] \left\{ \begin{array}{l} W^0, \\ W^{j+1}= W^j+ \Delta t \frac{f_u(t^{j},W^{j})+f_u(t^{j+1},W^{j+1})}{2}, \qquad j=0,1,...,J-1 \end{array} \right. [/math]

3. Runge-Kutta

[math] \left\{ \begin{array}{l} W^0, \\ W^{j+1}= W^j+ \frac{\Delta t}{6} \; (k_1+k_2+k_3+k_4), \qquad j=0,1,...,J-1 \\ k_1=f_u(t^{j},W^{j}) \\ k_2=f_u(t^{j}+\frac{\Delta t}{2},W^{j}+\frac{k_1}{2}\Delta t)\\ k_3=f_u(t^{j}+\frac{\Delta t}{2},W^{j}+\frac{k_2}{2} \Delta t)\\ k_4=f_u(t^{j}+\Delta t,W^{j}+k_3 \Delta t) \end{array} \right. [/math]

3 MATLAB code

%% Wave equation approximation: Method of lines-Euler in time
 clear all
 
%% Datos del problema
a=0; b=2; T=4;
 % Dato inicial y condiciones frontera
f=@(x,t) (x+t)*0;
u0=@(x) exp(-12*(x-1).^2);
u1=@(x) 24*(x-1).*exp(-12*(x-1).^2);
g1=@(t) 0*sin(t);
g2=@(t) 0*sin(t);

%% Discretización
N=20;   % number of space subintervals
mu=1; % relation between space and time discretization
h=(b-a)/N; % space step
xi=a+h:h:b-h; % interior space nodes (only interior nodes have equations with Dirichlet BC)
x=a:h:b;  % space nodes (including boundary ones)
dt=mu*h; % time step
t=0:dt:T;  % vector of times

%% Escribimos el sistema: U''=-K*U+F(t), U(0)=U0, U'(0)=U1. Calculamos K, U0, U1 y F
K=2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=K/h^2; % aproximation of -u_xx: matrix (N-1) x (N-1)
U0=u0(xi)'; % vector with the initial position
U1=u1(xi)'; % vector with the initial velocity
G1=zeros(N-1,1); G1(1)=1; % to include the boundary term in the first equation
G2=zeros(N-1,1); G2(N-1)=1; % % to include the boundary term in the last equation 
F=@(t) f(xi,t)'+g1(t)/h^2*G1+g2(t)/h^2*G2; % vector F

%% Escribimos el sistema: W'=-A*W+G(t), W(0)=W0. Calculamos A, W0 y G
A=[zeros(N-1),-eye(N-1);K,zeros(N-1)];
W0=[U0;U1];
G=@(t) [zeros(N-1,1);F(t)];

%% Métodos Explícitos: calculamos la función fu para escribir el sistema: W'=fu(t,W)
fu=@(t,W) -A*W+G(t); 

%% Bucle en tiempo
sol(1,:)=u0(x); % Vamos guardando la solución en la matriz sol (N+1)x(length(t))
W=W0; % Creamos un vector U que va a tener la aproximación en cada instante de tiempo. Inicialmente U=U0.
for k=1:length(t)-1
    % W=W+dt*fu(t(k),W); % Método de Euler en tiempo 
    % Pare resolver por el método del Trapecio descomentar la línea siguiente:
    W=(eye(2*N-2)+dt*A/2)\((eye(2*N-2)-dt*A/2)*W+dt*(G(t(k))+G(t(k+1)))/2);
    sol(k+1,:)=[g1(t(k+1)),W(1:N-1)',g2(t(k+1))]; % guardamos la solución en la fila correspondiente de la matriz sol.
end

%% Dibujamos
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)

4 Results

Consider the particular case: [math](a,b)=(0,2),\; u_0(x)=e^{-12(x-1)^2}, \; u_1(x)=24(x-1)e^{-12(x-1)^2}, \; g_1(t)=g_2(t)=0, \; f(x,t)=0, h=1/10 [/math]

left

--Carlos Castro (discusión) 15:09 31 ene 2013 (CET)