Całkowanie numeryczne f(x)dt – jednoczesne iterowanie dwóch wartości

0

Hej,
mam za zadanie napisać program, w którym zaimplementuję wzór:

W=f(x)+const*Integral(f2(V))dx

Chodzi tutaj o całkowanie np. metodą prostokątów. Mój problem polega na tym, że w funkcji wykonującej całkowanie nie wiem jak uwzględnić to, że x przyjmuje różne wartości (od -80, do 20), a całkujemy po t. Chciałabym, żeby program działał w pętli, gdzie iteracja byłaby po t, np. od 0 do 35.
Mój kod programu:

#include <iostream>
#include <math.h>

using namespace std;

double f2(double V){
    double vs=0;
    if(V<=-23){
        vs=0;
    }
    if(V>-23){
        vs= 2*pow(10,-7)*pow(V,6)-9*pow(10,-6)*pow(V,5)-4*pow(10,-4)*pow(V,4)+0.0088*pow(V,3)+0.1685*pow(V,2)-4.4818*V-64.826;
    }
    return vs;
}


double f(double x){
    double a=0;
    if(x<=5){
        a=0;
    }
    if(x>=5 && x<=10){
        a=(x-5)/5;
    }
    if(x>=10 && x<=15){
        a=1;
    }
    if(x>=15 && x<=20){
        a=(20-x)/5;
    }
    if(x>=20){
        a=0;
    }
    return a;
}

float calkowanie(float xp, float xk, int n){
    
    float calka;
    float h;

    h=(xp-xk)/(float)n;//krok
    
    calka=0;
    for(int i=1; i<=n; i++){
        calka+=fun2(xp+i*h)*h;
    }
    
    return calka;
    
}

int main(){
    
    
    double vv=0;
    double stala=7;
    
    for(int i=0; i<50; i++){
     
        vv=f(i)+stala*calkowanie(0,30,10);
            cout<<vv<<endl;
     
    return 0;
}

Proszę o pomoc - jak zmienić funkcję wykonującą całkowanie/ pętle for tak, aby iterowały się zarówno wartości V jak i x?

0

Ktoś już rozpisał w Kompendium całkowanie numeryczne metodą prostokątów i nawet podał przykład w kodzie ;)

Poza tym chyba troszkę nie rozumiem, podajesz w formułach f(x), a mówisz o całkowaniu i iterowaniu po t nie mówiąc, czym ma być to t - liczbą podprzedziałów, na który chcesz dzielić przedział całkowania?

EDIT

jak zmienić funkcję wykonującą całkowanie/ pętle for tak, aby iterowały się zarówno wartości V jak i x?

W funkcji calkowanie obliczasz wyrażenie xp + i*h, gdzie h=dlugosc_przedzialu/liczba_podprzedzialow, i przekazujesz to jako argument do fun2 (notabene definiujesz funkcję f2, wywołujesz fun2), gdzie pełni rolę v. Zatem v zmienia się gdy iterujesz w funkcji calkowanie po i. W czym więc problem?

0
superdurszlak napisał(a):

Ktoś już rozpisał w Kompendium całkowanie numeryczne metodą prostokątów i nawet podał przykład w kodzie ;)

Poza tym chyba troszkę nie rozumiem, podajesz w formułach f(x), a mówisz o całkowaniu i iterowaniu po t nie mówiąc, czym ma być to t - liczbą podprzedziałów, na który chcesz dzielić przedział całkowania?

EDIT

jak zmienić funkcję wykonującą całkowanie/ pętle for tak, aby iterowały się zarówno wartości V jak i x?

W funkcji calkowanie obliczasz wyrażenie xp + i*h, gdzie h=dlugosc_przedzialu/liczba_podprzedzialow, i przekazujesz to jako argument do fun2 (notabene definiujesz funkcję f2, wywołujesz fun2), gdzie pełni rolę v. Zatem v zmienia się gdy iterujesz w funkcji calkowanie po i. W czym więc problem?

ok, skoro xp+i*h odpowiada mojemu V, to w którym miejscu kodu całkuję i iteruję po x?

1
zozolek napisał(a):

Hej,
mam za zadanie napisać program, w którym zaimplementuję wzór:

W=f(x)+const*Integral(f2(V))dx

Nie rozumiem tego zapisu. FYI forum wspiera TEX-a: W(x) = f_1(x) + c \cdot\int_a^b f_2(x, y)dy
Nie wiem, czy o to chodzi, popraw tak, żeby było wiadomo jaki jest wzór.

zozolek napisał(a):

Chodzi tutaj o całkowanie np. metodą prostokątów. Mój problem polega na tym, że w funkcji wykonującej całkowanie nie wiem jak uwzględnić to, że x przyjmuje różne wartości (od -80, do 20), a całkujemy po t. Chciałabym, żeby program działał w pętli, gdzie iteracja byłaby po t, np. od 0 do 35.

Podaj wzór to może będzie jaśniejsze o co chodzi.

Samo całkowanie metodą prostokątów można zrobić tak:

template<typename T, typename F>
auto integrate(T a, T b, size_t n, F f) -> decltype(a * f(a))
{
    decltype(f(a)) sum {};
    auto dx = (b - a) / n;
    for (size_t i = 0; i < n; ++i)
        sum += f(a + 0.5 * dx + i * dx);

    return sum * dx;
}

auto result = integrate(0.0, 35.0, 1024, f);
0

W=f(x)+c*\int_a^b f2(V) dx
gdzie moje V ma przyjmować wartości od -80 do 20, a x od 0 do 35

kod programu:

#include <iostream>
#include <math.h>

using namespace std;

double f2(double V){
    double vs=0;
    if(V<=-23){
        vs=0;
    }
    if(V>-23){
        vs= 2*pow(10,-7)*pow(V,6)-9*pow(10,-6)*pow(V,5)-4*pow(10,-4)*pow(V,4)+0.0088*pow(V,3)+0.1685*pow(V,2)-4.4818*V-64.826;
    }
    return vs;
}


double f(double x){
    double a=0;
    if(x<=5){
        a=0;
    }
    if(x>=5 && x<=10){
        a=(x-5)/5;
    }
    if(x>=10 && x<=15){
        a=1;
    }
    if(x>=15 && x<=20){
        a=(20-x)/5;
    }
    if(x>=20){
        a=0;
    }
    return a;
}

float calkowanie(float xp, float xk, int n){
    
    float calka;
    float h;

    h=(xp-xk)/(float)n;//krok
    
    calka=0;
    for(int i=1; i<=n; i++){
        calka+=f2(xp+i*h)*h;
    }
    
    return calka;
    
}

int main(){

    double vv=0;
    double stala=7;
    
    for(int i=0; i<35; i++){
     
        vv=fun4(i)+stala*calkowanie();
            cout<<vv<<endl;
     }  
    return 0;
}

chcę, żeby jednocześnie w pętli zmieniało się V, a całkowało się po x (którego nie ma we wzorze tej funkcji - funkcja zależy tylko od V), stąd mój problem - jak w tym przypadku numerycznie całkować po x

0

Zmienną całkowania ma być x. A w kodzie właśnie nie wiem jak uwzględnić to, żeby zmieniało się V i jednocześnie całkowało po x (x nie występuje we wzorze tej funkcji, którą chcę całkować)

1

Możesz to napisać poprawnie od strony matematycznej? Obecnie można to przekształcić tak:
W(x) = f(x) + c\cdot(b - a)\cdot f_2(V)
ale wątpię, żeby o to chodziło, bo całka staje się zbędna.

A może jeszcze lepiej przepisz dokładną treść zadania jakie masz zrobić, bo twojej interpretacji problemu to raczej nie zrozumiem.

0

Dokładnie zadaniem jest rozwiązanie równania:

dx/dt = dy/dt + f(x)```

0
zozolek napisał(a):

Dokładnie zadaniem jest rozwiązanie równania:

dx/dt = dy/dt + f(x)```

A które spośród \frac{dy}{dt}, \frac{dx}{dt}, f(x(t)) są podane w zadaniu?

Jeżeli wzory na

\frac{dy}{dt}, f(x)

są podane, to możesz przyjąć początkowe t, x(t) i nie całkować numerycznie, bo tym się nie rozwiązuje równań różniczkowych, tylko zastosować jedną z metod numerycznego rozwiązywania równań różniczkowych, np. Eulera (najmniej dokładna, ale chyba najprostsza):

x(t + \Delta{}t) = x(t) +  \Delta{}t \cdot{} g(x, t) = x(t) +  \Delta{}t \cdot{} \frac{\delta{}x}{\delta{}t}

Czyli w każdym kroku liczysz (w przybliżeniu ofc) pochodną x po t i do x dodajesz zmianę x (liczoną jako pochodną pomnożoną przez długość kroku).

0

podane jest y(t) i f(x(t)), nie znamy dx/dt

0

no to bajka, mając y(t) możesz albo obliczyć dy/dt analitycznie, jeśli nie jest zbyt skomplikowane, albo numerycznie i wykorzystać do obliczania dx/dt, które wykorzystujesz do obliczania x(t + dt), liczysz kolejny krok etc... Euler albo R-K wyższego rzędu i cyk, dwójeczka :)

0

Ok, policzyłam dy/dt, jest to po prostu stała (0, 1/5, -1/5 w zależności od przedziału - wartości t). Napisałam program liczący metodą RK 4go rzędu, z tym, że zwraca on jedną wartość. Ja w tym przypadku chciałabym mieć to w funkcji czasu, tak żeby na wykresie wyszedł przebieg w czasie. W jaki sposób i w którym miejscu w kodzie mogłabym uwzględnić iterację, żeby zrobić to w pętli for?

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

using namespace std;

float dvdt(float t, float v){
    return 50+2.82*pow(10,6)*sqrt(8.39*pow(10,30))*1.602*pow(10,-19)*pow(10,10)*5*exp(0.534*pow(10,-19)*v);
}

float rungeKutta(float x0, float y0, float x, float h){
    int n=(int)((x-x0)/h);
    float k1, k2, k3, k4;
    float y=y0;
    for(int i=1; i<=n; i++){
        k1=h*dvdt(x0, y);
        k2=h*dvdt(x0+0.5*h, y+0.5*k1);
        k3=h*dvdt(x0+0.5*h, y+0.5*k2);
        k4=h*dvdt(x0+h, y+k3);
        //cout<<y<<endl;
        y=y+(1.0/6.0)*(k1+2*k2+2*k3+k4);
        
        x0=x0+h;
    }
    return y;
}

int main() {
    
    float x0=0;
    float y=1;
    float x=2;
    float h=0.001;
    
    cout<<rungeKutta(x0, y, x, h)<<endl;
    
    return 0;
}

0

Chyba najprościej byłoby wypisać wartość funkcji (i ew. czas) po każdej iteracji - drukowanie wartości zwracanej przez funkcję rungeKutta możesz pominąć. Pętlę już masz - w rungeKutta. Zresztą masz już nawet zakomentowane cout << y;, wystarczy odkomentować. Jeśli chcesz uwzględnić iterację, możesz drukować także indeks czy wartość czasu:

float rungeKutta(float x0, float y0, float x, float h){
    int n=(int)((x-x0)/h);
    float k1, k2, k3, k4;
    float y=y0;
    cout<<"x,y(x)"<<endl;
    cout<<x0<<","<<y<<endl;
    for(int i=1; i<=n; i++){
        k1=h*dvdt(x0, y);
        k2=h*dvdt(x0+0.5*h, y+0.5*k1);
        k3=h*dvdt(x0+0.5*h, y+0.5*k2);
        k4=h*dvdt(x0+h, y+k3);
        y=y+(1.0/6.0)*(k1+2*k2+2*k3+k4);
 
        x0=x0+h;
        cout<<x0<<","<<y<<endl;
    }
    return y;
}

Jak chcesz robić z tego wykres, może wygodniej będzie użyć std::ofstream by zapisać wyniki do pliku np. CSV zamiast pisać do std::cout.

0

gdy tak robiłam otrzymywałam zależność liniową, a z tego równania powinno mi wyjść coś takiego (calculated surface potential) Zrzut ekranu 2018-12-3 o 18.47.22.png

0

poprawiłam jeszcze kilka rzeczy i obecnie kod wygląda tak

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

using namespace std;

double E=8.854*pow(10,-12);
double e=1.602*pow(10,-19);
double m=9.109*pow(10,-31);
double n=1*pow(10,10);
double kbT=3;
double A=20;
double f=10;
double k=4;
double df=pow(10,-4);



float dydt(float t){
    if(t<=5){
        return 0;
    }
    if(t>5 && t<=10){
        return 1/5;
    }
    if(t>10 && t<=15){
        return 0;
    }
    if(t>15 && t<=20){
        return -1/5;
    }
    else{
        return 0;
    }
}
float dvdt(float t, float v){
    if(v<=-23){
        return dydt(t);
    }
    else{
        return dydt(t)+df/(k*E)*2*pow(10,-7)*pow(v,6)-9*pow(10,-6)*pow(v,5)-4*pow(10,-4)*pow(v,4)+0.0088*pow(v,3)+0.1685*pow(v,2)-4.4818*v-64.826;
    }
}
/*
float dvdt(float t, float v){
    return dydt(t)+df/(k*E)*e*sqrt(8*kbT/(3.14*m))*n*(A/4)*exp(e*v/kbT);
}
*/
float rungeKutta(float x0, float y0, float x, float h){
    int n=(int)((x-x0)/h);
    float k1, k2, k3, k4;
    float y=y0;
    for(int i=1; i<=n; i++){
        k1=h*dvdt(x0, y);
        k2=h*dvdt(x0+0.5*h, y+0.5*k1);
        k3=h*dvdt(x0+0.5*h, y+0.5*k2);
        k4=h*dvdt(x0+h, y+k3);
        cout<<y<<endl;
        y=y+(1.0/6.0)*(k1+2*k2+2*k3+k4);
        
        x0=x0+h;
    }
    return y;
}

int main() {
   
    
    float x0=0;
    float y=0; //poczatkowa wartosc w x0
    float x=2;
    float h=0.001;
    
    rungeKutta(x0, y, x, h);
    
    return 0;
}

natomiast wartości, które otrzymuję na wyjściu: output.txt
wykres nie jest już liniowy, ale to wciąż nie to co chciałam otrzymać.

1
zozolek napisał(a):

poprawiłam jeszcze kilka rzeczy i obecnie kod wygląda tak

float dydt(float t){

Staraj się unikać używania typu float - mamy 2018 rok, naprawdę możesz sobie pozwolić na zmarnowanie 8B zamiast 4B na pojedynczą liczbę zmiennoprzecinkową. float ma tragiczną precyzję, do bodajże 6-go miejsca po przecinku, lepiej używać typów podwójnej precyzji niż potem męczyć się z kosmicznymi błędami zaokrągleń :)

Zauważ, że z jednej strony operujesz na dość dużych liczbach, rzędu +/- 1e+1, z drugiej stosunkowo małymi, np. rzędu 1e-4. Praktycznie wyczerpujesz tym precyzję float. O wartościach rzędu 1e-12 czy 1e-19 to już nie wspomnę - nawet double ma z grubsza 16 cyfr po przecinku precyzji.

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