Soluciones Numéricas de la Ecuación del Telegrafista

Question

Actualmente estoy trabajando en un breve informe sobre la Ecuación del Telegrapher para mi curso de cálculo fraccionario. Todavía soy nuevo en las Ecuaciones del Telegrapher, pero sé que se utilizan para describir señales eléctricas que viajan a lo largo de un cable de transmisión (ya sea un cable coaxial, una microcinta, etc).

En resumen, derivé la Ecuación del Telegrapher al analizar los componentes elementales de una línea de transmisión:

\begin{eqnarray} \frac{\partial^2 v}{\partial t^2} + \frac{(LG+RC)}{LC}\frac{\partial v}{\partial t} + \frac{RG}{LC}v & = & \frac{1}{LC}\frac{\partial^2 v}{\partial x^2} \end{eqnarray}

Donde $v(x,t) = v$ es el voltaje a través del trozo de alambre, $R$ es la resistencia, $L$ es la inductancia, $G$ es la conductancia y $C$ es la capacitancia. Para simplificar, definiré $c= 1/\sqrt{LC}$, $a = c^2(LG+RC)$ y $b = c^2RG$, entonces obtenemos: \begin{eqnarray} \frac{\partial^2 v}{\partial t^2} + a\frac{\partial v}{\partial t} + bv & = &c^2\frac{\partial^2 v}{\partial x^2} \end{eqnarray}

Aquí es donde surge mi pregunta, todavía soy nuevo en la resolución numérica de EDPs, y mucho menos en la resolución de EDPs fraccionarias.

¿Podría alguien ayudarme a entender mejor cómo resolver este problema numéricamente (sin tener en cuenta las derivadas fraccionarias por ahora)? Supongo que necesitamos condiciones de contorno. Lamentablemente, no sabría cuáles serían las condiciones de contorno de una Ecuación del Telegrapher. Supongo que para mantener las cosas más simples, podemos asumir que estamos tratando con un sistema sin pérdidas, por lo tanto $R=G =0$. Así que nuestra EDP se simplifica mucho a:

\begin{eqnarray} \frac{\partial^2 v}{\partial t^2} & = &c^2\frac{\partial^2 v}{\partial x^2} \end{eqnarray}

que podemos usar para considerar la siguiente aproximación:

\begin{eqnarray} \frac{v(x,t+\Delta t)-2v(x,t)+v(x,t-\Delta t)}{\Delta x}& = & c^2\frac{v(x+\Delta x,t)-2v(x,t)+v(x-\Delta x,t)}{\Delta x} \end{eqnarray} Luego resolviendo para el tiempo futuro obtenemos:

\begin{eqnarray} v(x,t+\Delta t)& = & r^2[v(x+\Delta x,t)+v(x-\Delta x,t)] + 2(1-r^2)v(x,t)-v(x,t-\Delta t) \end{eqnarray}

Ahora tengo una idea aproximada de cómo configurar la solución en MATLAB, pero como soy nuevo en MATLAB y no tengo idea de cuál podría ser una condición de contorno de una ecuación del telegrapher (ya que no es mi campo de estudio), estoy terriblemente atascado y no he logrado avanzar tanto como me gustaría.

---

De todos modos, te agradezco sinceramente todo el tiempo que has tomado para leer esta publicación. También te agradezco de antemano por tu contribución, ayuda, retroalimentación y más.

Que tengas un día maravilloso.

Answer 1

No puedo ayudar con las derivadas fraccionarias, pero sí con las EDP "normales". Lo que suelo hacer es utilizar las herramientas de ODE de Matlab para encargarme de los pasos temporales y solo discretizar las derivadas espaciales. Entonces, utilizando $v_i(t)$ para denotar la solución en el $i$-ésimo punto de la rejilla espacial, obtienes ecuaciones como esta:

$$ \frac{\partial^2v_i}{\partial t^2}+a\frac{\partial v_i}{\partial t}+bv_i=\frac{c^2}{\Delta x^2}\left(v_{i-1}-2v_i+v_{i+1}\right). $$

Obviamente podrías utilizar un stencil diferente si lo deseas. También debes pensar en las condiciones de contorno, supondré que son $v_0=\alpha$ y $v_N=\beta$.

Esto te da un sistema de ecuaciones diferenciales matriciales. Olvidando el lado izquierdo por un momento, nota que podemos escribir el lado derecho como una multiplicación de matrices:

$$\frac{c^2}{\Delta x^2}\begin{bmatrix}\alpha & 0 & 0 & 0 & \ldots & 0 \\ 1 & -2 & 1 & 0 & \ldots & 0 \\ 0 & 1 & -2 & 1 & \ldots & 0 \\ 0 & 0 & 0 & 0 & \ldots & \beta\end{bmatrix}\begin{bmatrix}v_0\\v_1\\ \vdots\\v_N\end{bmatrix},$$ y a partir de ahora llamaré a la matriz de diferenciación $M$.

Ahora debes escribir esto como un sistema de EDO de primer orden, lo cual es bastante simple. Sea $u_{i,1}=v_i$ y $u_{i,2}=v'_i$, entonces $u'_{i,1}=u_{i,2}$ y

$$ u'_{i,2}=\frac{c^2}{\Delta x^2}(M\vec{u})_i-au_{i,2}-bu_{i,1}.$$ Escribir estas ecuaciones como vectores es más ordenado:

$$\vec{u}'_1=\vec{u}_2$$ y $$ \vec{u}'_2=\frac{c^2}{\Delta x^2}M\vec{u}_1-au_2-bu_1.$$

Puedes resolver esto utilizando ode45 en Matlab, con condiciones iniciales que tendrás que calcular. No es trivial, pero aquí está el código que resuelve $$\frac{\partial^2v}{\partial t^2}+a\frac{\partial v}{\partial t}+bv=c\frac{\partial^2v}{\partial x^2}. $$

a=1;b=1;c=1; %// coeficientes de la EDP
dx=0.05; %// espaciado de la rejilla espacial
x=(0:dx:1).'; %// rejilla espacial
N=length(x); 
%// matriz de diferenciación, utilizando diferencias centradas en puntos interiores
M=diag(ones(1,N-1),-1)+diag(-2*ones(1,N))+diag(ones(1,N-1),1);
M(1,1:3)=[1 -2 1]; %// y diferenciación hacia adelante/atrás en los límites
M(end,(N-2):N)=[1 -2 1];

%// Nuestro vector u es [u1;u2], por lo que la derivada tiene que ser [u2;M*u1-u2-u1]
%// más corrección para las condiciones de contorno
dudt=@(t,u)[[-c*sin(c*t);u(N+2:end-1);-c*sin(1+c*t)];
    c^2*M*u(1:N)/(dx^2)-a*u(N+1:end)-b*u(1:N)];

options=odeset('MaxStep',dx); %// algún tipo de condición CFL, no lo trabajé
                              %// correctamente, pero esto funciona así que está bien!

[t,u]=ode45(dudt,[0 15],[cos(x);-c*sin(x)],options); %// pasos temporales con ode45
                        %// la condición inicial es [v;dv/dt] en t=0

for i=1:length(t) %// Para ver los resultados
    plot(x,u(i,1:N))
    axis([0 1 -1 1])
    drawnow
    pause(0.01)
end

Si estableces $a=b=0$ entonces puedes comparar la solución numérica con la solución exacta $\cos(x+ct)$.

¡Espero que eso ayude!