Algorytm inercji I rzędu

Question

Witam, mam problem z napisaniem algorytmu inercji I rzędu.

Wzór jest następujący:

Mój algorytm na podstawie tego wzoru:

double Object::calculate(double value)
{
  y1 = y2;
  y2 = out;
  out = (mK * value) - mT * (y1 - y2);
  return out;
}

Konstruktor klasy Object przyjmuje 2 argumenty K i T

Object::Object(double K, double T)
{
   mK = K;
   mT = T;
   y1 = 0;
   y2 = 0;
   out = 0;
}

Pętla z maina do sprawdzania wartości:

    Object tempObject = Object(1, 1);
    double val = 20;
    for(int i = 0; i < 30; i++)
    {
        std::cout << i << ". " << "temp = " << tempObject.value() << std::endl;
        tempObject.calculate(val);
    }

Wartości z konsoli:

0. temp = 0
1. temp = 20
2. temp = 40
3. temp = 40
4. temp = 20
5. temp = 0
6. temp = 0
7. temp = 20
8. temp = 40
9. temp = 40
10. temp = 20
11. temp = 0
12. temp = 0
13. temp = 20
14. temp = 40
15. temp = 40
16. temp = 20
17. temp = 0
18. temp = 0
19. temp = 20
20. temp = 40
21. temp = 40
22. temp = 20
23. temp = 0
24. temp = 0
25. temp = 20
26. temp = 40
27. temp = 40
28. temp = 20
29. temp = 0

Teoretyczne (tak jak powinien działać) sprawdzenie tego członu w MatLab:

Jak widać pomimo wzmocnienia K = 1 i stałej czasowej T = 1 wartości oscylują wokół tej zadanej. Wydaje mi się, że problem leży w tej całce, ew. stałej czasowej. Ale nie mam na to pomysłu. Może ktoś miał podobny problem ?

Comments

Dlaczego założyłeś, że dt==1? Takie równania należy liczyć tak, bo kolejne wyniki było minimalnie różne. Jak widać MathLab zanim doszedł do 20 obliczył całkiem sporo wyników.

---

tak w klasie, maja wartości 0

Answer 1

Nie będę dociekał, jak to powinno wyglądać w C++. Tym bardziej nie będę się wdawał w klasy (można mnie nazwać przeciwnikiem OOP, ale używanie klas w momencie, kiedy wystarczy jedna prosta funkcja, to już przykład overengineering).

Natomiast jeśli chodzi o obliczenie wyjścia na podstawie wejścia, to mam następujące uwagi:

0) W Twoich wzorach jest założenie, że warunek początkowy jest zerowy:

y(0) = 0

1) Dalej jest potrzebna jakaś funkcja do całkowania numerycznego.

Bo mamy przecież:

dy(t) / dt = K/T * [ x(t) - y(t) ]

U Ciebie jest jakiś inny wzór. W sumie to nie wiem, skąd się wziął. (A może autor miał coś zupełnie innego na myśli?)

No to dla najprostrzego całkowania, metodą Eulera, może być tak (kod w Pythonie, można sprawdzić działanie w JDoodle, jak ktoś nie ma zainstalowanego):

def fPrim(x, K, T):
  def prim(y):
    return K/T * (x - y)
  
  return prim

def euler(f, start = 0, end = 1, h = 0.125, initial = 0):
  y = initial
  t = start
  out = [y]
  while t < end:
    y = y + h * f(y)
    out.append(y)
    t = t + h
  return out

print( euler(fPrim(20, 1, 1), end=5) )

Funkcja euler przyjmuje jako pierwszy parametr funkcję zwracającą pochodną po czasie dla danej wartości y.

Poniżej coś podobnego w C (C++ się tykał nie będę ;-)), testowane w JDoodle:

#include<stdio.h>

double deriv(double x, double K, double T, double y) {
  return K/T * (x - y);
}

int main() {
  double t = 0;
  double y = 0;
  double h = 0.125;
  while ( t < 5 ) {
    printf( "%.3f %.2f\n", t, y );
    y = y + h * deriv(20, 1, 1, y);
    t = t + h;
  }
  printf( "%.3f %.2f\n", t, y );

  return 0;
}

Zwracam przy tym drobną uwagę na krok całkowania: h = 0.125. Można zastosować inny, ale dla bezpieczeństwa (numerycznego) powinna być to wartość 1 / 2^N (N - całkowite). Ja zastosowałem N = 3, czyli 1/8.

Przy całkowaniu numerycznym można też zastosować inne metody, np. Rungego-Kutty - będą zapewnie dawały mniejsze błędy.