[Fortran] Całkowanie numeryczne

Yassamet
Użytkownik
Użytkownik
Posty: 36
Rejestracja: 29 cze 2013, o 14:16
Płeć: Kobieta
Lokalizacja: Wrocław
Podziękował: 3 razy
Pomógł: 1 raz

[Fortran] Całkowanie numeryczne

Post autor: Yassamet »

Witajcie.
Mam taką całkę: \(\displaystyle{ \int_{1}^{\infty} t^{-n} e^{-xt} dt = E_{n} (x)}\)
Jest to funkcja, której zadaje się z góry parametry \(\displaystyle{ n}\) i \(\displaystyle{ x}\), a ona zwraca wartość. Potrzebuję ją zaprogramować numerycznie, zdecydowałam się na metodę trapezów wedle wzoru:
\(\displaystyle{ \frac{h}{2} \sum_{i=1}^{ \infty } (f_{x_{i}} + f_{x_{i+1}} )}\)
Gdzie \(\displaystyle{ h}\) oznacza krok między zmiennymi, a funkcje \(\displaystyle{ f}\) oznaczają funkcję podcałkową. Zaprogramowałam i właśnie - liczy, ale źle. Stąd moje pytanie, czy to pomysł jest zły, czy może w kodzie jest błąd? Metoda działała dobrze dla prostych całek typu \(\displaystyle{ x^{-2}}\).
Załączam także kod i pozdrawiam serdecznie.

Kod: Zaznacz cały

subroutine Expo(n,x,E)
      real :: E, A, B1, B2, x, t, zmienna
      integer :: n, i
      t = 1.0
      E = 0.0
      i = 0
      A = (t**(-n))*(exp(-x*t))
      zmienna = 0.0001 * A
      do
      i = i+1
      B1 = (t**(-n))*(exp(-x*t))  
      t = t + 0.01
      B2 = (t**(-n))*(exp(-x*t))
      E = E + B1 + B2
        if(B2.lt.zmienna) then
        exit
        end if
      end do
      write(*,*) i
      E = E * 0.01 / 2
      end subroutine
      
Ostatnio zmieniony 8 cze 2016, o 07:56 przez Afish, łącznie zmieniany 1 raz.
Powód: Poprawa wiadomości.
SlotaWoj
Użytkownik
Użytkownik
Posty: 4211
Rejestracja: 25 maja 2012, o 21:33
Płeć: Mężczyzna
Lokalizacja: Kraków PL
Podziękował: 2 razy
Pomógł: 758 razy

[Fortran] Całkowanie numeryczne

Post autor: SlotaWoj »

Yassamet pisze:Zaprogramowałam i właśnie - liczy, ale źle.
A co to znaczy źle?
Zrealizowałem Twój algorytm w Excelu (bo tak było najprościej) i otrzymałem \(\displaystyle{ E_2(\pi)=8,951464\cdot10^{-3}}\) z błędem względnym \(\displaystyle{ 7,825\cdot10^{-5}}\) w stosunku do wartości dokładnej.

Mam następujące uwagi: procedura jest mało optymalna i nie podoba mi się warunek zakończenia pętli (\(\displaystyle{ f(t)/f(1)<0,0001}\) – dlaczego właśnie taki?).
Yassamet
Użytkownik
Użytkownik
Posty: 36
Rejestracja: 29 cze 2013, o 14:16
Płeć: Kobieta
Lokalizacja: Wrocław
Podziękował: 3 razy
Pomógł: 1 raz

[Fortran] Całkowanie numeryczne

Post autor: Yassamet »

Ja porównywałam wynik z programu z portalu wolphramalfa, to w wielu przypadkach różnica była rzędu wyniku lub nawet wyższa.
Jeśli chodzi o ten warunek, to po prostu przyjmuję taką dokładność. Funkcja jest w całej dziedzinie zbieżna, zatem nie obawiam się dużego błędu z tego tytułu.

Dodam też, że napisałam tę procedurę raz jeszcze, poczynając od realizacji dla prostszych funkcji i już generuje dobre wyniki. Prawdopodobnie gdzieś powyżej jest literówka.
SlotaWoj
Użytkownik
Użytkownik
Posty: 4211
Rejestracja: 25 maja 2012, o 21:33
Płeć: Mężczyzna
Lokalizacja: Kraków PL
Podziękował: 2 razy
Pomógł: 758 razy

[Fortran] Całkowanie numeryczne

Post autor: SlotaWoj »

  1. A co to znaczy w wielu wypadkach? Jakie \(\displaystyle{ n}\) i \(\displaystyle{ x}\) dawało złe wyniki?
  2. Przepisywałem ww. kod na Pascal i żadnej literówki nie zauważyłem.
  3. Wróć do pierwotnego kodu (tego, który daje złe wyniki) i prześledź go (np. wypisując i analizując wartości pośrednie).
  4. Przetestuj wpływ na wynik deklaracji typów real/double precision.
  5. Jeśli poprawny kod daje błędne wyniki, to albo kompilator Fortranu jest „do bani”, albo z komputerem coś się dzieje.
ODPOWIEDZ