Komentarze (0)  Wizyty (1 623)  Aktualnie osób (3)  Tutaj (1)

User: Hasło:

 Numerical Methods of Inverse Laplace Transform

Reklamy Google

Istnieją różne numeryczne metody obliczania odwrotnej transformaty Laplace'a, takie jak np :

Metoda oparta na całkowaniu równania odwrotnej transformaty - bardzo wolna ale dokładna

Algorytm Zakian'a - bardzo prosta w użyciu, szybka w implementacji ale nie nadająca się do układów oscylacyjnych

Metody oparta na FFT - szybka w obliczeniach i dokładna ale wymagająca algorytmu FFT

Metoda z użyciem szeregów Laguerre'a

Algorytm Walsa

W oparciu o pozostałości w biegunach Postaram się przedstawić kilka powyższych implementacji i podać przykłady.

Metoda oparta na całkowaniu równania transformaty odwrotnej

Metoda ta jest bardzo prosta w implementacji - wprost z definicji - i od niej też zacznę. Polega ona na obliczaniu całek dla każdej chwili czasu t wg. wzoru :

Poniżej przedstawiam przykładowa funckję dla Delphi obliczającą dla chwili czasu t wartość odwrotnej transformaty. W przykładzie uprościłam tak metodę całkowania żeby byla ona prosta w analizie. Użyłam tutaj najprostszej metody całkowania - prostokątów :

function Integration(t, Omega, Sigma, nCalk : Double) : Double;
var Delta, wi, fi : Double;
    val           : TComplex;
    i             : Integer;
begin
    wi := 0;
    Suma := 0;
    Delta := Omega / nCalk;

    for i:=0 to Round (nCalk) do
    begin
        wi := wi + Delta;
        val := RealToCplx (0,wi*t);

        f2 := MulCplx (ExpCplx(val),K(RealToCplx(Sigma,wi))).Re;
        Suma := Suma + Delta * fi;
    end;

    Result := Suma * EXP (Sigma * t) / pi;
end;

Omega jest maksymalną częstotliwością wyrażoną w [ rad / sek ]. Ja generalnie uzywam wartości Omega = 400 - wystarcza to aby także wysokie częstotliwości weszły w skład wykresu. Sigma to mała wartość dodatnia w zakresie około 0.01 - 1. nCalk to liczba kroków całkowania dla każdej chwili czasu. Ja ustawiam ok. 400.

W funkcji Integrate użyta jest inna funkcja - K ( Z : TComplex ) : TComplex - zwraca ona wartość w postaci zespolonej transmitancji K ( s ). Do obliczeń wartości funkcji pobieranej z łańcucha znaków używam własnej roboty parsera matematycznego dla liczb zespolonych - znacznie to ułatwia sprawę.

Odpowiedź impulsowa ekstrapolatora zerowego rzędu obliczona za pomocą metody całkowania równania odwrotnej transformaty -
widać oscylacje przy narastającym zboczu.

---

Algorytm Zakian'a

Metoda ta jest najprostszą sposród dostępnych do zaimplementowania. Polega ona na na obliczaniu sum ważonych dla każdej interesującej nas chwili czasu t w ten sposób:

Jak widać do obliczeń potrzebne są nam tylko stałe Alpha ( i ) oraz K ( i ). Oznaczenia:

f ( t ) - poszukiwana funkcja z dziedziny czasu

F ( s ) - dana funkcja z dziedziny operatora Laplace'a (oryginał Laplace'a)

Alpha ( i ) - stałe brane z popniższej tablicy

K ( i ) - stałe brane z popniższej tablicy

t - parametr czasu

REAL( x ) - część rzeczywista liczby x

Tabela zespolonych współczynników Alpha ( i ) oraz K ( i ) dla kolejnych pięciu i :

i               Alpha ( i )                       K ( i )
--------------------------------------------------------------------
1      12.83767675 + j 1.666063445      -36902.08210 + j 196990.4257
2      12.22613209 + j 5.012718792      +61277.02524 - j 95408.62551
3      10.93430308 + j 8.409673116      -28916.56288 + j 18169.18531
4      8.776434715 + j 11.92185389      +4655.361138 - j 1.901528642
5      5.225453361 + j 15.72952905      -118.7414011 - j 141.3036911

Jak widac metoda jest prosta i szybka obliczeniowo. Jednak należy pamiętać o tym, że nie można za jej pomocą wyznaczyć wartości f ( t ) dla chwili czasu t = 0. Warto również zaznaczyć że metoda ta nie nadaje się dla układów oscylacyjnych.

Proszę zobaczyć przebiegi uzyskane za pomoca Zakian'a oraz komentarze do nich :

Układ oscylacyjny 2-go rzędu. Wykres biały (metoda Zakiana) - widać że po 5 punkcie przegięcia metoda już nie działa dobrze.

Inercja 1-go rzędu. Wykres biały - Zakian - widać poprawne działanie metody.

Visitors

Reklamy Google

---

Komentarze (0)

Zapraszam do pisania komentarzy w tym temacie.

Napisz komentarz

Pamiętaj, aby wypełnić wszystkie pola.

Temat: Podpis: