Ver el código de «Diferencias finitas para la ecuación de ondas»


== Planteamiento: ecuación de ondas ==

Consideramos el siguiente problema

<math>\begin{equation}\label{diferencia1}\left\{  \begin{array}{lll} u_{tt}-u_{xx}=\frac{4e^{-t}}{1+x^2}, \hspace{0.4cm} x \in (0,1), \;\; t>0, \\ u_x(0,t)=2te^{-t^2}, \;\; u(1,t)=\frac{t}{1+t^2}, \hspace{0.4cm} t>0, \\ u(x,0)=x^2(1-x), \;\; u_t(x,0)=x, \hspace{0.4cm} x \in(0,1). \end{array} \right.\end{equation}
</math>
Vamos a aproximar la solución por diferencias finitas en <math> x\in (0,1) </math>  y <math> t \in (0,3) </math>.


== Discretización espacial ==

Primeramente introducimos una partición del intervalo <math> x\in (0,1) </math>. Definimos el número de subintervalos N=5 y consideramos la malla en espacio $x_i=a+idx, \; i=0,1,2,3,4,5$. Llamaremos $dx=\frac{b-a}{N}=\frac{1}{5}$ la distancia entre dos puntos consecutivos de la malla.
En los puntos <math> x_i </math> tenemos

<math> u_{tt}(x_i,t)-u_{xx}(x_i,t)=\frac{4e^{-t}}{1+x_i^2}, \hspace{0.4cm} i=1,2,3,4. </math>

Reemplazamos <math>-u_{xx}(x_i,t)</math> por <math>\frac{-u(x_{i-1},t)+2u(x_i,t)-u(x_{i+1}t)}{dx^2},</math>
obteniendo entonces una aproximación de la ecuación en derivadas parciales
<math> \begin{equation}\label{diferencia2}u_{tt}(x_i,t)+\frac{-u(x_{i-1,t}+2u(x_i,t)-u(x_{i+1}t)}{dx^2}=\frac{4e^{-t}}{1+x_i^2}, \hspace{0.4cm} i=1,2,3,4.\end{equation} </math>

Si llamamos $u_i(t)=u(x_i,t), \; i=0,1,2,3,4,5$ el sistema (\ref{diferencia2}) lo podemos escribir en la forma
<math>\begin{equation}\label{diferencia3}\left\{ \begin{array}{llll} u''_1(t)+\frac{-u_0(t)+2u_1(t)-u_2(t)}{dx^2}=\frac{4}{1+x_1^2}e^{-t}  \\  u''_2(t)+\frac{-u_1(t)+2u_2(t)-u_3(t)}{dx^2}=\frac{4}{1+x_2^2}e^{-t}  \\  u''_3(t)+\frac{-u_2(t)+2u_3(t)-u_4(t)}{dx^2}=\frac{4}{1+x_3^2}e^{-t}  \\  u''_4(t)+\frac{-u_3(t)+2u_4(t)-u_5(t)}{dx^2}=\frac{4}{1+x_1^2}e^{-t}  \end{array}  \right.\end{equation}
</math>

No conocemos el valor de $u(0,t)$, solo el valor de la derivada con respecto a $x$,

<math>u_x(0,t)= \lim_{\Delta x \rightarrow 0} \frac{u(0+\Delta x,t)-u(0,t)}{\Delta x}=2te^{-t^2}.</math> 

Tomamos $\Delta x= dx$ y aproximamos $\lim_{\Delta x \rightarrow 0} \frac{u(0+\Delta x,t)-u(0,t)}{\Delta x}$ por $\frac{u(dx,t)-u(0,t)}{d x}$. Entonces

<math>2te^{-t^2}=\frac{u(dx,t)-u(0,t)}{d x}=\frac{u(x_1,t)-u(0,t)}{d x}=\frac{u_1(t)-u_0(t)}{d x} \Rightarrow u_0(t)= u1(t)-2te^{-t^2}dx,</math>

y el sistema \ref{diferencia3} podemos escribirlo en la forma
<math>\begin{equation}\label{diferencia4b}\left\{ \begin{array}{llll} u''_1(t)+\frac{-u_1(t)+2te^{-t^2}dx+2u_1(t)-u_2(t)}{dx^2}=\frac{4}{1+x_1^2}e^{-t}  \\  u''_2(t)+\frac{-u_1(t)+2u_2(t)-u_3(t)}{dx^2}=\frac{4}{1+x_2^2}e^{-t}  \\  u''_3(t)+\frac{-u_2(t)+2u_3(t)-u_4(t)}{dx^2}=\frac{4}{1+x_3^2}e^{-t}  \\  u''_4(t)+\frac{-u_3(t)+2u_4(t)-\frac{t}{1+t^2}}{dx^2}=\frac{4}{1+x_1^2}e^{-t}  \end{array}  \right.\end{equation}</math>
o en forma equivalente
<math>\begin{equation}\label{diferencia4}\left\{ \begin{array}{llll} u''_1(t)+\frac{-u_1(t)+2u_1(t)-u_2(t)}{dx^2}=\frac{4}{1+x_1^2}e^{-t} -\frac{2te^{-t^2}}{dx} \\  u''_2(t)+\frac{-u_1(t)+2u_2(t)-u_3(t)}{dx^2}=\frac{4}{1+x_2^2}e^{-t}  \\  u''_3(t)+\frac{-u_2(t)+2u_3(t)-u_4(t)}{dx^2}=\frac{4}{1+x_3^2}e^{-t}  \\  u''_4(t)+\frac{-u_3(t)+2u_4(t)}{dx^2}=\frac{4}{1+x_1^2}e^{-t}+\frac{t}{1+t^2} \frac{1}{dx^2} \end{array}  \right.\end{equation}</math>

Si llamamos 

<math> U(t)= \left(   \begin{array}{llll} u_1(t) \\ u_2(t) \\ u_3(t)  \\ u_4(t) \end{array}  \right), \qquad  K= \left(  \begin{array}{cccc} 1 & -1 & 0 & 0 \\-1 & 2 & -1 & 0 \\    0 & -1 & 2 & -1 \\   0 & 0 & -1 & 2 \; \end{array} \right), \qquad  F(t)=\left( \begin{array}{llll}  \frac{4}{1+x_1^2}e^{-t} -\frac{2te^{-t^2}}{dx} \\   \frac{4}{1+x_2^2}e^{-t} \\ \frac{4}{1+x_3^2}e^{-t} \\  \frac{4}{1+x_4^2}e^{-t}  +\frac{t}{1+t^2} \frac{1}{dx^2}  \end{array}  \right), </math>

el sistema \ref{diferencia4}, junto con la condici\'on inicial, puede escribirse en la forma

<math>\begin{equation}\label{diferencia5}\left\{  \begin{array}{ll}  U''(t)=-KU(t)+F(t)    \\ U(0)=\left( \begin{array}{llll} x_1^2(1-x_1) \\ x_2^2(1-x_2) \\ x_3^2(1-x_3) \\ x_4^2(1-x_4) \end{array}      \right) , \;\; U'(0)=\left( \begin{array}{llll} x_1 \\ x_2 \\ x_3 \\ x_4 \end{array}      \right)  \end{array}  \right.    .\end{equation}</math>


== Discretización temporal ==

El siguente paso es resolver este sistema de 4 ecuaciones diferenciales de segundo orden.Empezamos primeramente utilizando el m\'etodo de Euler. Esto implica pasar de este sistema de 4 ecuaciones diferenciales de segundo orden a un sistema de 8 ecuaciones diferenciales de primer orden . Introducimos las funciones $u'_1(t)=v_1(t), \; u'_2(t)=v_2(t), \; u'_3(t)=v_3(t), \; u'_4(t)=v_(t),$, las nuevas icognitas son $8 $ funciones que las designamos por $u_1(t)$, $u_2(t)$, $u_3(t)$, $u_4(t)$, $v_1(t)$, $v_2(t)$, $v_3(t)$, $v_4(t)$ que satisfacen el sistema de 8 ecuaciones diferenciales de primer orden
<math>\begin{equation}\label{diferencia6}\left\{ \begin{array}{llllllll} u_1'(t)=v_1(t) \\ u_2'(t)=v_2(t) \\ u_3'(t)=v_3(t) \\ u_4'(t)=v_4(t) \\  v'_1(t)=-\frac{-u_1(t)+2u_1(t)-u_2(t)}{dx^2}+\frac{4}{1+x_1^2}e^{-t} -\frac{2te^{-t^2}}{dx} \\  v'_2(t)=-\frac{-u_1(t)+2u_2(t)-u_3(t)}{dx^2}+\frac{4}{1+x_2^2}e^{-t}  \\  v'_3(t)=-\frac{-u_2(t)+2u_3(t)-u_4(t)}{dx^2}+\frac{4}{1+x_3^2}e^{-t}  \\  v'_4(t)=\frac{-u_3(t)+2u_4(t)}{dx^2}+\frac{4}{1+x_1^2}e^{-t}+\frac{t}{1+t^2} \frac{1}{dx^2}  \end{array} \right.\end{equation}
</math>
LLamando $V(t)= \left(   \begin{array}{llll} v_1(t) \\ v_2(t) \\ v_3(t)  \\ v_4(t) \end{array}  \right)$, \ref{diferencia5} es equivalente a
<math>\begin{equation}\label{diferencia7}\left\{ \begin{array}{ll}U'(t)=V(t)  \\V'(t)=-KU(t)+F(t)  \\ U(0)=\left( \begin{array}{llll} x_1^2(1-x_1) \\ x_2^2(1-x_2) \\ x_3^2(1-x_3) \\ x_4^2(1-x_4) \end{array}      \right) , \;\; V(0)=\left( \begin{array}{llll} x_1 \\ x_2 \\ x_3 \\ x_4 \end{array}      \right)  \end{array}  \right.\end{equation}</math>
Si tomamos un paso temporal $dt$, definimos $t_j=0+jdt$, $U_n$ y $V_n$ designa la soluci\'on aproximada de $U(t_n)$ y $V(t_n)$ respectivamente,   el m\'etodo de Euler aplicado a \ref{diferencia7} queda en la forma
<math>\left(  \begin{array}{ll} U_{n+1}  \\ V_{n+1}  \end{array} \right) =\left(  \begin{array}{ll} U_{n}  \\ V_{n}  \end{array} \right)+dt \left(  \begin{array}{ll} U'_{n}(t_n)  \\ V'_{n}(t_n)  \end{array} \right)=\left(  \begin{array}{ll} U_{n}  \\ V_{n}  \end{array} \right)+dt\left(  \begin{array}{ll} V_{n}  \\ -KU_n+F(t)  \end{array} \right)</math> 


o escrito en otra forma
<math>\begin{equation}\label{diferencia8}\left\{ \begin{array}{lll} U_{n+1}=U_n+dtV_n   \\ V_{n+1} = V_n+dt(-KU_n+F(t_n))    \\ U_0=\left( \begin{array}{llll} x_1^2(1-x_1) \\ x_2^2(1-x_2) \\ x_3^2(1-x_3) \\ x_4^2(1-x_4) \end{array}      \right) , \;\; V_0=\left( \begin{array}{llll} x_1 \\ x_2 \\ x_3 \\ x_4 \end{array}      \right)  \end{array}  \right.\end{equation}</math>

== Programa para OCTAVE o MATLAB ==

{{matlab|codigo=
%%%% Diferencias finitas para resolver la ecuacion de ondas
clear all
L1=0; L2=1;       % extremos del intervalo
%% Discretizacion espacial
N=10;             % numero de subintervalos 
dx=(L2-L1)/N;     % paso en espacio
x=L1+dx:dx:L2-dx;
U=(x.^2.*(1-x))'; % posicion inicial
V=x';             % velocidad inicial
TF=3;             % tiempo final
%% discretización en tiempo
dt=dx/10;         % paso en tiempo
t=0:dt:TF;
%% matriz que aproxima -u_xx
K=(1/(dx^2))*(2*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1));
K(1,1)=1/dx^2;
g=4./(1+x.^2)';       % segundo miembro
sol(1,:)=[U(1),U',0]; % matriz que contiene la solucion para dibujarla
f=zeros(N-1,1);       % vector con la aportacion de la frontera 
for k=1:length(t)-1
    f(1)=-(2*t(k)*exp(-t(k)^2)/dx);  % condicion neumann en x=L1
    f(N-1)=(t(k)/((1+t(k)^2)*dx^2)); % condicion dirichlet en x=L2
    U=U+dt*V;                         % método explícito
    V=V+dt*(-K*U+g*exp(-t(k))+f);
    sol(k+1,:)=[U(1)-2*t(k)*exp(-t(k))*dx , U', t(k)/(1+t(k)^2)];
end
%Dibujo de la solucion
x1=L1:dx:L2
[X,T]=meshgrid(x1,t)
surf(X,T,sol)
}}

== Ejemplo ==

[[Archivo:Diferenciasfiinitaseuler.jpg|600px|thumb|centro|Aproximación numérica de la ecuación de ondas]]

[[Categoría:Grado en Ingeniería Civil y Territorial]]
[[Categoría:Ecuaciones Diferenciales]]