Laplace Inversion utilizando el teorema de residuos

Question

ACTUALIZACIÓN: he intentado responder a esta pregunta

Estoy tratando de invertir la transformada de Laplace usando Cauchy del Teorema de los Residuos. La transformación es

$$\Theta(s)=\frac{F(s)}{G(s)}$$

donde

$$ F(s)=[1+\gamma P(s)] [As^{-1}+\theta(0)] $$

y $$ G(s)=[1+\gamma P(s)] (s+A) + \gamma BP(s) $$

en que $\theta(0)$ es la condición inicial de la inversa de la función de transformación en el dominio del tiempo, $\gamma$, $A$, e $B$ son reales constantes, y $P(s)$ está dada por la serie

$$ P(s)=s\sum_{k=1}^{k=N} \frac{C_k}{s+\lambda_k^2} $$

donde $N \ge 50$. Los términos de $C_k$ e $\lambda_k$ son los verdaderos Fourier-Bessel coeficientes reales y valores propios, que son el resultado de otro problema acoplado.

Aquí está el método que he utilizado:

Hay un polo en $s=0$ y el restante $N$ polos $z_k, k=1\cdots N$ son los ceros de $G(s)$, todos de los cuales se encuentran en el Eje Real negativo, y cada uno de los cuales he calculado con residuos muy cerca de la máquina de redondeo. Siguiendo el Teorema de los Residuos, la transformada Inversa de Laplace puede ser escrito como

$$ \theta(t)=1+\sum_{k=1}^{k=N}e^{z_k t}F(z_k) \lim_{s\a z_k}\frac{s-z_k}{G(s)} $$

con lo cual el límite indeterminado puede ser resuelto con la Regla de L'Hôpital rendimiento

$$ \theta(t)=1+\sum_{k=1}^{k=N} e^{z_k t}\frac{F(z_k)}{G'(s)} $$

donde $G'(s)$ es la derivada de w.r.t. $s$ de $G(s)$

$$ G'(s)=[1+\gamma P(s)] +\gamma (s+a+B)P'(s) $$

y $P'(s)$ es la derivada de w.r.t. $s$ de $P(s)$

$$ P'(s)=\sum_{k=1}^{k=N} \frac{C_k}{s+\lambda_k^2}-s\sum_{k=1}^{k=N} \frac{C_k}{(s+\lambda_k^2)^2} $$

Temas:

He observado que la evaluación de la $\theta(t)$ , no parece muy justo, especialmente para los pequeños de veces y, en particular, para grandes valores de $\gamma$. Además, se demuestra una sensibilidad a los valores propios y los Coeficientes de $C_k$ que no es consistente con las expectativas de las observaciones. En lugar de inundar al lector con más detalles en esta etapa, estaré encantado de actualizar a la pregunta con la información más relevante siguientes comentarios.

Mi preguntas específicas:

Es el enfoque que yo he usado con el Teorema de los Residuos fundamentalmente correcto? ¿Hay algo más que debo ser consciente de la hora de resolver problemas del tipo de arriba?

Muchas gracias, como siempre.

Answer 1

Tal vez la incorporación de un pre-procesamiento para reducir la función de transferencia a su expresión más simple, hacer que el proceso de la inversión más precisa. Para $N = 2$, el final de la función de transferencia se presenta, adecuada para aplicar el método de

$$ \Theta(s) = \frac{A \left(C_1 \gamma s \left(\lambda _2^2+s\right)+\left(\lambda _1^2+s\right) \left(C_2 \gamma s+\lambda _2^2+s\right)\right)}{s \left(C_1 \gamma s \left(\lambda _2^2+s\right) (a+B+s)+\left(\lambda _1^2+s\right) \left(C_2 \gamma s (a+B+s)+(A+s) \left(\lambda _2^2+s\right)\right)\right)} $$

with denominator coeficients

$$ \left\{ \begin{array}{c} 0 \\ A \lambda _1^2 \lambda _2^2 \\ \left(\lambda _2^2+A+(A+B) \gamma C_2\right) \lambda _1^2+\left(A+(A+B) \gamma C_1\right) \lambda _2^2 \\ \lambda _1^2+\lambda _2^2+A+\gamma C_2 \left(\lambda _1^2+A+B\right)+\gamma C_1 \left(\lambda _2^2+A+B\right) \\ \gamma \left(C_1+C_2\right)+1 \\ \end{array} \right. $$

Follow a symbolic processing script in MATHEMATICA to implement this simplification

P[s_, n_] := s Sum[Subscript[C, k]/(s + Subscript[lambda, k]^2), {k, 1, n}]
F[s_, n_] := (1 + gamma P[s, n]) A/s
G[s_, n_] := (1 + gamma P[s, n]) (s + A) + gamma B P[s, n]
Phi[s_, n_] := F[s, n]/G[s, n]

n = 2;
Phisn = Phi[s, n] // Together
num = Phisn // Numerator;
den = Phisn // Denominator;
Phis0 = num/den // FullSimplify
den0 = Phis0 // Denominator

Attached a python script to obtain the inverse Laplace transform from a transfer function $G(s) = \frac{N(s)}{D(s)}$ with the $N(s) = a_n s^n + \cdots + a_0$ and $D(s) = b_m s^m+\cdots b_0$ The $N$ coefficients are given in num and the $D$ coefficients in den both in power decreasing style. $G(s)$ se supone estrictamente correcto. El guión da una secuencia de respuesta de los términos que deben ser añadidos a compuesto el momento de la inversión.

import scipy.signal as R 


num = [6.92462, 73.8839, 400.845, 1455.68, 3946.03, 8455.99, \
       14824.9, 21736.1, 27038.6, 28803.0, 26421.5, 20924.8, \
       14306.8, 8423.79, 4249.23, 1821.19, 655.104, 194.266, \
       46.253, 8.49133, 1.12494, 0.0951142, 0.00380817]
den = [1, 10.9544, 61.0185, 227.587, 634.139, 1398.5, 2527.38,\
       3827.46, 4929.55, 5452.66, 5212.0, 4320.11, 3108.84, 1940.74,\
       1048.1, 487.221, 193.437, 64.822, 18.0145, 4.04279, 0.702988, \
       0.0884995, 0.00711991, 0.000271338, 0]

#num = [3, 2, 2]
#den = [5, 2, 2, 1]

(r,p,k) = R.residue(num,den,tol=1e-5)

toler1 = 1e-5
toler2 = 1e-15
i      = 0
j      = 1

while i < len(r):
    alphai = real(r[i])
    betai  = imag(r[i])
    gammai = real(p[i])
    deltai = imag(p[i])

    if abs(betai) > toler1:
        if abs(alphai) < toler2:
            alphai = 0   
        print("{4} >>> 2*exp({0}*t)*({1}*cos({2}*t) + {3}*sin({2}*t))".format(gammai,alphai,deltai,-betai,j))
        i += 2
    else:
        if abs(gammai) > toler2:
            print("{2} >>> {0}*exp({1}*t)".format(1*alphai,1*gammai,j))
        else:
            print("{1} >>> {0}".format(1*alphai,j))
        i += 1
    j += 1

La inversión del problema da como resultado

1 >>> 14.0347831855
2 >>> -0.483039070616*exp(-0.187218948148*t)
3 >>> 2*exp(-0.309980077984*t)*(-1.97460596883*cos(-0.0673684262979*t)+0.0408608048522*sin(-0.0673684262979*t))
4 >>> 2*exp(-0.317487772418*t)*(-0.853151294864*cos(-0.193564461065*t)+0.11474944314*sin(-0.193564461065*t))
5 >>> 2*exp(-0.262710062176*t)*(-0.192942085124*cos(-0.380624549122*t)+-0.196136519663*sin(-0.380624549122*t))
6 >>> 2*exp(-0.0566781486718*t)*(-2.04346564851*cos(-0.473314360884*t)+-0.784055678847*sin(-0.473314360884*t))
7 >>> 2*exp(-0.253402873944*t)*(-0.224610305483*cos(-0.602941037526*t)+0.573864761592*sin(-0.602941037526*t))
8 >>> 2*exp(-0.18001778095*t)*(-0.214405624466*cos(-0.953393232567*t)+-0.105678203189*sin(-0.953393232567*t))
9 >>> 2*exp(-0.25943110475*t)*(-0.377299755441*cos(-1.02160634165*t)+-0.266212072126*sin(-1.02160634165*t))
10 >>> 2*exp(-1.07768853278*t)*(0.0995875511961*cos(-0.120122202754*t)+-0.456124154697*sin(-0.120122202754*t))
11 >>> 2*exp(-0.674142284507*t)*(-0.207698716256*cos(-1.01903055486*t)+-0.0654651577739*sin(-1.01903055486*t))
12 >>> 2*exp(-1.08508395079*t)*(-0.222022363057*cos(-0.651792664191*t)+-0.740528774445*sin(-0.651792664191*t))
13 >>> 2*exp(-0.906967936955*t)*(-0.565257846658*cos(-0.938474102998*t)+-0.35832504633*sin(-0.938474102998*t))

Que debe ser añadido para dar la respuesta a un impulso. En la figura se muestran dos curvas. Uno de los asumido exacta y la otra, la calculada por la secuencia de comandos.

Answer 2

He descubierto la causa de la discrepancia. La solución del problema subyacente en el dominio del tiempo es

$$ \theta(t)=1+\sum_{k=1}^{k=N} e^{z_k t}\frac{F(z_k)}{G'(s)} $$

como se ha indicado anteriormente. El problema es que en el momento $0$, la condición inicial $\theta(0))$ debe ser honrado. Esto requiere que $$ \theta(0)=1+\sum_{k=1}^{k=N} \frac{F(z_k)}{G'(s)} $$

que es violado por una pequeña cantidad de dinero cuyo preciso de la magnitud depende del tamaño de los valores de los primeros polos son. He encontrado que si he añadido esta diferencia a la solución, todo funciona a la perfección. Vea la figura de abajo. La Gaver-Stehfest enfoque utiliza una función algoritmo de muestreo y es muy adecuado para la difusión de los problemas.

Lo que todavía no consigue, aunque la verdad es que he incluido $\theta(0)$ en la expresión de $F(s)$, esto no es satisfecho en lo que se creo una integral definida de la evaluación en el Método de los residuos. ¿Por qué?!