Różniczkowanie metodą Strliinga

0

Witam i zgłaszam się do was jak co weekend z kolejnym działaniem. Tym razem na warsztat dostałem wzór Stirlinga. Dokładni opis mozna znaleźć tutaj: http://www.eti.pg.gda.pl/katedry/ksa/dydaktyka/Algorytmy_obliczeniowe//Wykl_AlgorOblicz_8.pdf

Więc mam programik, który wymyśliłem w C kiedyś i wygląda on tak:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

int i, j;
double avg[8], di[8][8], h, s;
double xi[8] = { 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4 }; // wspolczynniki x
double yi[8] = { 5.000000, 2.500000, 1.666667, 1.250000, 1.000000, 0.833333, 0.714286 }; // wspolczynniki y

int wyznacz_avg() {
    j = 0;
    for( i=0; i<sizeof( yi )/sizeof( yi[0] ); i++ ) {
        if( i % 2 != 0 ) {
            avg[j] = ( di[i][sizeof( yi )/sizeof( yi[0] )/2-i/2-1]+di[i][sizeof( yi )/sizeof( yi[0] )/2-i/2] )/2;
            j += 1;
        }
    }

    return( 1 );
}

int wyznacz_d() {
    for( i=0; i<sizeof( yi )/sizeof( yi[0] ); i++ ) {
        di[0][i] = yi[i];
    }

    for( j=0; j<sizeof( yi )/sizeof( yi[0] )-1; j++ ) {
        for( i=0; i<sizeof( yi )/sizeof( yi[0] )-j-1; i++ ) {
            di[j+1][i] = di[j][i+1] - di[j][i];
        }
    }

    return( 1 );
}

int wyznacz_h() {
    h = xi[1] - xi[0];

    return( 1 );
}

int stirling() {
    s = ( 1/h ) * ( avg[0] - avg[1]/6 + avg[2]/30 );

    return( 1 );
}

int main() {
    if( wyznacz_h() && wyznacz_d() ) {
        if( wyznacz_avg() ) {
            if( stirling() ) {
                printf( "wynik rownania: %f\n", s );
            }
            else {
                printf( "blad w trakcie obliczania rownania\n" );
            }
        }
        else {
            printf( "blad w trakcie srednich\n" );
        }
    }
    else {
        printf( "blad w trakcie obliczania wartosci kroku oraz delty\n" );
    }

    return( 0 );
}
 

Teraz mam to napisac w pythonie, wiec nie wiem czy dobrze zrobiłem ale zacząłem to jakoś przerabiać. na razie jestem tutaj:

'''
Created on 16-11-2013

@author: Konrad
'''

#!/usr/bin/python

xi = [ 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4 ]; 
yi = [ 5.000000, 2.500000, 1.666667, 1.250000, 1.000000, 0.833333, 0.714286 ]; 
avg = [];
h = 0.0;

def wyznacz_sred():
    for xi,yi in zip(xi,yi):
    

i kompletnie nie wiem jak sobie poradzic z tymi tablicami dwuwymiarowymi. Ogólnie liczenie h można sobie podarować i traktować to jako stałą, muszę mieć średnią, d no i potem można podstawiać do wzoru tak jak tutaj w funkcji stirling. Proszę o pomoc.

1

Wzór Stirlinga kojarzy się raczej z szacowaniem silnii, ale biorąc to co linkowałeś.

Mam problem ze zorientowaniem się co liczy Twój kod w C :P.
Jeśli chcesz wyliczyć pochodną to najpierw przybliżasz funkcję f z tej tabeli przez jakąś funkcję (różniczkowalną) F (w tym przypadku wielomian stirlinga), po czym ją różniczkować i liczyć że pochodna F dobrze przybliża pochodną f.

Noo i chodzi o to że Twój kod wypisuje jedną wartość (wynik rownania: -1.056549), więc na pewno cała pochodna to nie jest ;]. Najwyżej pochodna w jakimś punkcie, ale jeśli tak to nie mogę się zorientować w jakim.

(wg. wzorów, F(x) = ... i dF(x)/dx = ... - czyli dla każdego x mamy jakąś wartość (a nie jedną wartość dla całej funkcji))

Więc może zapytam co chcesz w zasadzie policzyć? :P

0

Więc tak: powiedzmy mam zbiór punktów:

xi[8] = { 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4 }; // wspolczynniki x
yi[8] = { 5.000000, 2.500000, 1.666667, 1.250000, 1.000000, 0.833333, 0.714286 }; // wspolczynniki y

i chcę to policzyć w ten sposób:

Najpierw otrzymać taką tablicę róźnic centralnych(tablica przykładowa 1 z brzegu) :

user image

Potem podstawiam do wzoru:

user image

gdzie h - to krok czyli x[1] - x[0]
a reszta to te wspolczynniki z tablicy

Jak widać bo tam nie dopisalem x i y mamy podane a kolumne delta f(x) obliczamy odejmując ostani element od przedostaniego czyli np: 0,18232 - 0,09531 i tak az do konca.

1

Najpierw otrzymać taką tablicę róźnic centralnych(tablica przykładowa 1 z brzegu) :

Na pewno o tym samym PDF (http://www.eti.pg.gda.pl/katedry/ksa/dydaktyka/Algorytmy_obliczeniowe//Wykl_AlgorOblicz_8.pdf) mówimy? :P Takiej tabeli tam nigdzie np. nie ma.

Ale pomijając to -

http://i.imgur.com/avuIbyT.jpg

To już bardziej pasuje do tego co Twój kod robi - tylko że

  1. s = ( 1/h ) * ( avg[0] - avg[1]/6 + avg[2]/30 ); we wzorze są trzy kropki na końcu, to może (i ma) mieć więcej niż trzy wyrazy
msm napisał(a)

Noo i chodzi o to że Twój kod wypisuje jedną wartość (wynik rownania: -1.056549), więc na pewno cała pochodna to nie jest

wzór: mówi f'(x) = coś tzn. pochodna funkcji opisywanej przez tabelę w punkcie x jest równa coś. Możesz wyliczyć tą wartość w jakimś punkcie najwyżej.

U Ciebie:
avg[j] = ( di[i][sizeof( yi )/sizeof( yi[0] )/2-i/2-1]+di[i][sizeof( yi )/sizeof( yi[0] )/2-i/2] )/2;
I niezbyt to się do wzoru ma (przede wszystkim nie ma nigdzie x, więc nie może to byc ta sama funkcja).

0

No nie z tego pdf'a. Ten co podałem to bardzie ogolne jest opisane a tutaj podalem juz to co odnosi sie do samego programu.

  1. Spójrz że we wzorze bierzemy same "wpolczynniki nieparzyste". Czyli z tej tabelki bierzemy fi oraz fi3. Poniewaz to tablica roznic cenrtalnych to wybieramy "wspolczynnik ze srodka" Tutaj mamy 4 wyrazy(liczba parzysta) wiec musimy je do siebie dodac. Dzialamy:

0,095310 + 0,10535 / 2 = 0,10033

tak samo dla fi3:

0,00175 + 0,00237 / 2= 0,00206 i podstawiamy do wzoru:

user image

to 0.1 w tym wzorze bierze sie z tego ze h czyli krok wynosi 0.1. Obliczamy to z tabelki z X - 1 - 0.9 = 0.1

  1. To własnie ta funkcja sprawdza i jezeli liczba jest parzysta to robi to co opisalem wyzej.
1
  1. to pisałem odnośnie kodu, chyba że już założyłeś tam że tablice zawsze są takiej samej wielkości (ale wtedy z kolei niepotrzebne te sizeof( yi )/sizeof( yi[0] )).

  2. Rozumiem wszystko, tylko nie rozumiem dlaczego. Masz funkcję:

user image
f(x) = 1/h * ( (df(x - 1/2h) + df(x + 1/2h)) / 2 - ... )
czyli wartość funkcji dla argumentu x jest równa 1/h * ( (fd(x - 1/2h) + df(x + 1/2h) ) / 2 ... ).
Tzn. dlaczego bierzesz wartość ze środka tabeli, zamiast z okolic x jak chce ten wzór?

0

Tam wyżej ci tłumaczyłem. Dokładnie tez się nie znam na tym ale na zajeciach tak więc powtarzam i to jest na bank dobrze. Ta tabela roznic centralnych która ci tam wrzuciłem wyżej wygląda jak wygląda. We wzorze jak widzisz bierzemy wartości dla fi nieparzystego i podstawiamy do tego f(x) = 1/h(fi1 + 1/6* fi3) bo tak wynika z tabelki. Jakby w tabeli w wyniku tych odejmowań wyszło jeszcze fi5 wtedy wzór byłby inny. Teraz mamy tylko fi oraz fi3.

Wiec po naszym liczeniu mamy

f(x) = 1/h(0,10033 + 1/6*0,00206) i jest wynik. Jeżeli w tych tabelkach:

user image
user image

mielibysmy nieparzysta ilosc wspolczynnikow np. 3 to zawsze bralibysmy ten srodkowy. Mam nadzieje ze wytłumaczyłem ci to o co Ci chodziło.

2

Tam wyżej ci tłumaczyłem. Dokładnie tez się nie znam na tym ale na zajeciach tak więc powtarzam i to jest na bank dobrze

Właśnie ja wyżej wytłumaczyłem dlaczego nie jest dobrze, a przynajmniej nijak ma się do wzoru/materiałów które podałeś.

Cóż, zasadniczo jeśli chcesz to nie ma problemu w bezmyślnym przepisaniu tego (przetestowane tylko dla tego jednego przypadku, ale powinno być ok):

from math import floor

xi = [ 0.2,      0.4,      0.6,      0.8,      1.0,      1.2,      1.4      ]; 
yi = [ 5.000000, 2.500000, 1.666667, 1.250000, 1.000000, 0.833333, 0.714286 ]; 

def diffs(tab):
    for i in range(len(tab) - 1):
        yield tab[i + 1] - tab[i]

def diff_table(tab):
    while len(tab) > 1:
        tab = list(diffs(tab))
        yield tab

def get_avg(table):
    for i in range(0, len(table), 2):
        mid, mid2 = int((len(table[i]) - 1) / 2), int(len(table[i]) / 2)
        avg = (table[i][mid] + table[i][mid2]) / 2
        yield avg

def stirling():
    h = xi[1] - xi[0]
    table = list(diff_table(yi))
    avg = list(get_avg(table))
    return (1.0 / h) * (avg[0] - avg[1] / 6 + avg[2] / 30)

print stirling()

Uwaga 1)

def diffs(tab):
    for i in range(len(tab) - 1):
        yield tab[i + 1] - tab[i]

działa w uproszczeniu podobnie jak:

def diffs(tab):
    result = []
    for i in range(len(tab) - 1):
        result.append(tab[i + 1] - tab[i])
    return result

tylko jest ładniejsze i zwraca generator zamiast listy (czyli dalej zamiast tab = diffs(tab) piszę tab = list(diffs(tab))).
Jeśli przeszkadza Ci ten zapis to możesz zmienić wszystkie yieldy (pierwszy) na ręczne budowanie listy (drugie).

  1. Wynik tego kodu jest inny niż Twojego w C++, dlatego że Twój kod w C++ źle liczy indeksy do avg (przesunięte o jeden w prawo).

1 użytkowników online, w tym zalogowanych: 0, gości: 1