GNUPLOT aus einer Programmiersprache heraus nutzen

Antworten
asterix
Beiträge: 241
Registriert: 23.05.2018, 08:24

16.08.2019, 14:34

Auf der Suche nach einer brauchbaren Library zur Diagrammdarstellung in c++ bin ich auf Gnuplot gestossen. Das Programm lässt sich ziemlich einfach in beliebige Programmiersprachen einbinden. Die Kommunikation erfolgt über eine Pipe. Wenn man nur Daten nach Gnuplot schaufeln will, ist die Methode äußerst simpel zu realisieren:

Code: Alles auswählen

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main(int argc, char **argv)
  {
  FILE *gp;
  if ((gp = popen("gnuplot -persist","w"))==NULL)  /* pipe descriptor */
    exit(0);
  fprintf(gp, "set style data lines\n");
  fprintf(gp, "set grid\n");
  fprintf(gp, "plot '-'\n");  // Daten  Start
  for (double x=0; x < 6.28; x+=6.28/100)
    fprintf(gp, "%g %g  with lp linestyle 1\n", x,sin(x));  // Koordinatenpaar an Gnuplot senden
  fprintf(gp, "e\n");       // Daten Ende
  fflush(gp);
  fprintf(gp,"exit\n");
  fclose(gp);
  return 0;
  }

Das kleine c-Programm startet Gnuplot und überträgt die Daten die in einem xy-Diagramm dargestellt werden sollen. Die Skalierung erfolgt voreinstellungsmäßig automatisch. Die Ausgabe sieht dann so aus:

Abb1.gif
Abb1.gif (18.87 KiB) 184 mal betrachtet


Dass das Ganze dann noch mit dem altbekannten Algorithmus getestet wurde, versteht sich von selbst.

Abb2.gif
Abb2.gif (33.88 KiB) 184 mal betrachtet

Mit dem TDMGCC 5.1.0 kommt man mit der gleichen Konfiguration wie beim Java/Python-Test auf eine Laufzeit von etwa einer Sekunde.
Das ist meine persönliche Meinung dazu. Basierend auf einer nach bestem Wissen und Gewissen recherchierten Faktenlage.
asterix
Beiträge: 241
Registriert: 23.05.2018, 08:24

23.08.2019, 13:23

Das Quellpprogamm für die untere Grafik hat noch gefehlt.

Code: Alles auswählen

include <cstdlib>
#include <iostream>
#include <math.h>
#include <time.h>
#define PI 3.1415926535897932

class Vektor
  {
  public:
  double x;
  double y;
  double z;
  Vektor()
    {
    x=0;
    y=0;
    z=0;
    }
  Vektor(double x,double y,double z)
    {
    this -> x=x;
    this -> y=y;
    this -> z=z;
    }
/*
  Vektor(Vektor v)
    {
    this -> x=v.x;
    this -> y=v.y;
    this -> z=v.z;
    }
*/
  double Punktprodukt(Vektor v)
    {
    return this -> x * v.x + this -> y * v.y + this -> z * v.z ;
    }
  Vektor Kreuzprodukt(Vektor v)
    {
    return Vektor((this -> y * v.z) - (this -> z * v.y),(this -> z * v.x) - (this -> x * v.z),(this -> x * v.y) - (this -> y * v.x));
    }
  Vektor Set(Vektor v)
    {
    this -> x = v.x;
    this -> y = v.y;
    this -> z = v.z;
    return Vektor(this -> x,this -> y,this -> z);
    }
  Vektor Rot_x(double phi)
    {
    double y=this -> y*cos(phi)-this -> z*sin(phi);
    double z=this -> z*cos(phi)+this -> y*sin(phi);
    return Vektor(this -> x,y,z);
    }
  Vektor Rot_y(double phi)
    {
    double x=this -> x*cos(phi)-this -> z*sin(phi);
    double z=this -> z*cos(phi)+this -> x*sin(phi);
    return Vektor(x,this -> y,z);
    }
  Vektor Rot_z(double phi)
    {
    double x=this -> x*cos(phi)-this -> y*sin(phi);
    double y=this -> y*cos(phi)+this -> x*sin(phi);
    return Vektor(x,y,this -> z);
    }
  void Rot_xIs(double phi)
    {
    double y=this -> y*cos(phi)-this -> z*sin(phi);
    double z=this -> z*cos(phi)+this -> y*sin(phi);
    this -> y=y;
    this -> z=z;
    }
  void Rot_yIs(double phi)
    {
    double x=this -> x*cos(phi)-this -> z*sin(phi);
    double z=this -> z*cos(phi)+this -> x*sin(phi);
    this -> x=x;
    this -> z=z;
    }
  void Rot_zIs(double phi)
    {
    double x=this -> x*cos(phi)-this -> y*sin(phi);
    double y=this -> y*cos(phi)+this -> x*sin(phi);
    this -> x=x;
    this -> y=y;
    }
#include <math.h>

  Vektor Addieren(Vektor v)
    {
    return Vektor(this -> x + v.x,this -> y + v.y,this -> z + v.z);
    }
  Vektor Subtrahieren(Vektor v)
    {
    return Vektor(this -> x - v.x,this -> y - v.y,this -> z - v.z);
    }
  Vektor Skalarmultiplikation(double s)
    {
    return Vektor(this -> x*s,this -> y*s,this -> z*s);
    }
  Vektor Skalardivision(double s)
    {
    return Vektor(this -> x/s,this -> y/s,this -> z/s);
    }
  double Betrag()
    {
    return sqrt(this -> x*this -> x+this -> y*this -> y+this -> z*this -> z);
    }
  Vektor normieren()
    {
    double l=Betrag();
    return Vektor(this -> x/l,this -> y/l,this -> z/l);
    }
  };

class Simulation
  {
  public:
   double muo ;
   int steps;
   int stepsR;
   double r1;
   double r2;
   double Spalt;
   double phi0;

   double I1;
   double I2;
   double R1;
   double R2;

   double *Mphi;
   double *Mphi1;
   double *Mphi2;
   double *IMphi;

  int initg;

Simulation(double r2,double r1,double R2,double R1,double Spalt,double phi0,double I2,double I1,int steps,int stepsR)
  {
  initg=0;
  muo = 4.0*PI*1.0e-7;

  this -> steps = steps; //integration steps Rotor
  this -> stepsR = stepsR;  // integration steps Shleife
  this -> r1=r1;// [mm]    21.0     outer radius rotor #1
  this -> r2=r2;// [mm]    30.0     outer radius rotor #2
  this -> Spalt=Spalt;// [mm]    1.0     gap between rotors
  this -> phi0=phi0;// [deg]    ?5°       angle shift between rotors

  this -> I1=I1;// [kA]     0.88     current in loop #1
  this -> I2=I2;// [kA]  0.88   current in loop #2
  this -> R1=R1;// [mm]     5.0     radius loop #1
  this -> R2=R2;// [mm]     5.0     radius loop #2

  Mphi=NULL;
  Mphi1=NULL;
  Mphi2=NULL;
  IMphi=NULL;
  }


Simulation()
  {
  initg=0;
  muo = 4.0*PI*1.0e-7;
  steps = 512; //integration steps Rotor
  stepsR = 32;  // integration steps Shleife
  r1=0.021;// [mm]    21.0     outer radius rotor #1
  r2=0.030;// [mm]    30.0     outer radius rotor #2
  Spalt=0.004;// [mm]    1.0     gap between rotors
  phi0=-5;// [deg]    ?5°       angle shift between rotors

  I1=880;// [kA]     0.88     current in loop #1
  I2=880;// [kA]  0.88   current in loop #2
  R1=0.005;// [mm]     5.0     radius loop #1
  R2=0.005;// [mm]     5.0     radius loop #2

  Mphi=NULL;
  Mphi1=NULL;
  Mphi2=NULL;
  IMphi=NULL;

  }
~Simulation()
  {
  if(Mphi) delete(Mphi);
  if(Mphi1) delete(Mphi1);
  if(Mphi2) delete(Mphi2);
  if(IMphi) delete(IMphi);
  }

void init()
  {
  if(Mphi) delete(Mphi);
  if(Mphi1) delete(Mphi1);
  if(Mphi2) delete(Mphi2);
  if(IMphi) delete(IMphi);

  Mphi=new double[steps];
  Mphi1=new double[steps];
  Mphi2=new double[steps];
  IMphi=new double[steps];
  initg=1;
  }


Vektor Flussdichte_B1(Vektor Pb1,double phi2) //Biot Savart Flussdichte B1 verursacht durch S2
  {

  Vektor H(0,0,0);
  double dphidl2=2.0*PI/stepsR;
  Vektor Pdl2(r2,R2,0); //(1) S2 linke Schleife
  Pdl2=Pdl2.Rot_y(phi2);
  int n=1;
  while(n <= stepsR)
    {
    double phidl2=n*dphidl2;    // links herum
    Vektor Pdl2_nplus1(r2,R2*cos(phidl2),R2*sin(phidl2)); //(1) S2 linke Schleife
    Pdl2_nplus1=Pdl2_nplus1.Rot_y(phi2);
    Vektor dl2=Pdl2_nplus1.Subtrahieren(Pdl2); //(2)
    Vektor rd1=Pdl2.Subtrahieren(Pb1);  //(12)
    double abst=rd1.Betrag();
    Vektor dH=dl2.Kreuzprodukt(rd1).Skalardivision(abst*abst*abst);//()
    H=H.Addieren(dH);
    Pdl2=Pdl2.Set(Pdl2_nplus1);
    n+=1;
    }



  return H.Skalarmultiplikation(I2*muo/(4.0*PI));  //(5)
  }

Vektor Flussdichte_B2(Vektor Pb2,double phi1) //Biot Savart Flussdichte B2 verursacht durch S1
  {

  Vektor H(0,0,0);
  double dphidl1=2.0*PI/stepsR;
  Vektor Pdl1(-r1,R1,0); // (3)   S1 rechte Schleife
  Pdl1=Pdl1.Rot_y(phi1);
  int n=1;
  while(n <= stepsR)
    {
    double phidl1=n*dphidl1;    // links herum
    Vektor Pdl1_nplus1(-r1,R1*cos(phidl1),R1*sin(phidl1)); //(3)   S1 rechte Schleife
    Pdl1_nplus1=Pdl1_nplus1.Rot_y(phi1);
    Vektor dl1=Pdl1_nplus1.Subtrahieren(Pdl1); //(4)
    Vektor rd2=Pb2.Addieren(Pdl1);   //(6)
    double abst=rd2.Betrag();
    Vektor dH=dl1.Kreuzprodukt(rd2).Skalardivision(abst*abst*abst);//(5)
    H=H.Addieren(dH);
    Pdl1=Pdl1.Set(Pdl1_nplus1);
    n+=1;
    }
  return H.Skalarmultiplikation(I1*muo/(4.0*PI));  //(5)
  }

Vektor dphi_rotieren_S1(double phi2,double phi1,int n_) // Rotor-Teildrehung dphi
  {
     Mphi1[n_]=0; //(14)
     Vektor Fs(0,0,0);
     double dphidl1=2.0*PI/stepsR;
     Vektor Pdl1(-r1,R1,0); // (3)   S1 rechte Schleife
     Pdl1=Pdl1.Rot_y(phi1); // ()
     int n=1;
     while(n <= stepsR)
       {
       double phidl1=n*dphidl1;      // links herum
       Vektor Pdl1_nplus1(-r1,R1*cos(phidl1),R1*sin(phidl1)); //(3)   S1 rechte Schleife
       Pdl1_nplus1=Pdl1_nplus1.Rot_y(phi1);  //um Phi rotieren
       Vektor dl1=Pdl1_nplus1.Subtrahieren(Pdl1); //()
       Vektor Pb1(r1+r2+Spalt,0,0);
       Pb1=Pb1.Addieren(Pdl1); // (12)
       Vektor B1=Flussdichte_B1(Pb1,phi2); // (11)
       Vektor dFs1=dl1.Kreuzprodukt(B1).Skalarmultiplikation(I1); //(13)
       Fs=Fs.Addieren(dFs1);
       Mphi1[n_]+=(Pdl1.Kreuzprodukt(dFs1)).y; //(14)
       Pdl1=Pdl1.Set(Pdl1_nplus1);
       n+=1;
       }
  return Fs;
  }


Vektor dphi_rotieren_S2(double phi2,double phi1,int n_) // Rotor-Teildrehung dphi
  {
     Mphi2[n_]=0; //(14)
     Vektor Fs(0,0,0);
     double dphidl2=2.0*PI/stepsR;

     Vektor Pdl2(r2,R2,0); //(1) S2 linke Schleife
     Pdl2=Pdl2.Rot_y(phi2); // ()
     int n=1;
     while(n <= stepsR)
       {
       double phidl2=n*dphidl2;      // links herum
       Vektor Pdl2_nplus1(r2,R2*cos(phidl2),R2*sin(phidl2)); //(1) S2 linke Schleife
       Pdl2_nplus1=Pdl2_nplus1.Rot_y(phi2);  //um Phi rotieren
       Vektor dl2=Pdl2_nplus1.Subtrahieren(Pdl2); //()
       Vektor Pb2(r1+r2+Spalt,0,0);
       Pb2=Pb2.Subtrahieren(Pdl2); // (6)
       Vektor B2=Flussdichte_B2(Pb2,phi1); // (5)
       Vektor dFs2=dl2.Kreuzprodukt(B2).Skalarmultiplikation(I2); //(9)
       Fs=Fs.Addieren(dFs2);
       Mphi2[n_]+=(Pdl2.Kreuzprodukt(dFs2)).y; //(10)
       Pdl2=Pdl2.Set(Pdl2_nplus1);
       n+=1;
       }
  return Fs;
  }
void rotieren() // eine Umdrehung des Rotors
   {
   init();
   double dphi=2.0*PI/(double)steps;
   Vektor Fs(0,0,0);
   int n=0;
   while(n < steps)
     {
     double phi2=n*dphi;
     double phi1=-n*dphi+phi0/360.0*2.0*PI;
     phi2+=PI;
     phi1+=PI;
     Fs=dphi_rotieren_S2(phi2,phi1,n);
     Fs=dphi_rotieren_S1(phi2,phi1,n);
     Mphi[n]=Mphi1[n]-Mphi2[n];
   //  std::cout << n << "   " << Mphi[n]<< "   " <<  Mphi1[n]<< "   " <<  Mphi2[n] << "\n";
     ++n;
     }
   }
 };




FILE *plotpipe;
int open();
void plot(const char *x)
  {
  fprintf(plotpipe,"%s\n",x);
  }



int main(int argc, char *argv[])
  {
  int  steps = 512;
  int stepsR = 64;
  double r1=0.021;
  double r2=0.030;
  double Spalt=0.004;
  double phi0=-5;

  double I1=880;
  double I2=880;
  double R1=0.005;
  double R2=0.005;

  char param1[200];
  char param2[200];
  sprintf(param1,"\"steps=%d stepsR=%d R2=%.1fmm R1=%.1fmm  r2=%.1fmm r1=%.1fmm \"",steps,stepsR,R2*1000.0,R1*1000.0,r2*1000.0,r1*1000.0);
  sprintf(param2,"\"Spalt=%.1fmm Phi0=%.1f° I1=%.1fA I2=%.1fA \"",Spalt*1000.0,phi0,I1,I2);
  std::cout << param1 << "\n";
  std::cout << param2 << "\n";

  Simulation s(r2,r1,R2,R1,Spalt,phi0,I2,I1,steps,stepsR);
  clock_t  t0=clock();
  s.rotieren();
  clock_t  t1=clock();
  std::cout << "CLOCKS_PER_SEC=" << CLOCKS_PER_SEC << "   " << (t1-t0) << "  " << (t1-t0)*1000/CLOCKS_PER_SEC << "ms \n";

  if(!(plotpipe = popen("gnuplot -persist","w")))
    {
    printf("Error opening pipe to GNU plot. Check if you have it! \n");
    exit(0);
    }

  plot("set terminal wxt size 800,800 enhanced font 'Verdana,10' persist");
  plot("set multiplot layout 2,1");
  plot("set style line 1 lt rgb \"red\" lw 2");
  plot("set style line 2 lt rgb \"orange\" lw 2");
  plot("set object 1 rect from graph 0, graph 0 to graph 1, graph 1 back");
  plot("set object 1 rect fc rgb \"cyan\" behind");
  plot("set xlabel \"Winkel/°\"");
  plot("set style data  lines");
  plot("set grid");
  plot("set xrange [-180:180]");
  plot("set yrange [-0.009:0.0075]");

  plot("set title \"Drehmoment\"");
  plot("set label 8 at -170, -0.0068");
  plot("set label 9 at -170, -0.008");

  char param1_1[200];
  char param2_1[200];
  sprintf(param1_1,"set label 8 %s  tc lt 1",param1);
  plot(param1_1);
  sprintf(param2_1,"set label 9 %s  tc lt 1",param2);
  plot(param2_1);

  plot("set ylabel \"Drehmoment/Nm\"");
  plot("set key left reverse Left");

  plot("plot '-'  title \"M\" ls 1 ");       // plot '-' using ... with ...
  double dphi=2.0*PI/(double)steps;
  int n=0;
  while(n < steps)
    {
    char buf[100];
    double phi2=n*dphi;
    phi2-=PI;
    double x=phi2/(2.0*PI)*360;
    sprintf(buf,"%g %g",x, (float)s.Mphi[n]);
    plot(buf);
    ++n;
    }
  plot( "e");
  fflush(plotpipe);

  plot("set yrange [-0.25:0.03]");
  plot("set label 8 at -170, -0.21");
  plot("set label 9 at -170, -0.23");
  plot("set title \"Summe Drehmomente\"");
  plot("set ylabel \"Summe Drehmomente/Nm°\"");
  plot("plot '-'  title \"Int(M*dPhi)\" ls 2 ");       // plot '-' using ... with ...
  dphi=2.0*PI/(double)steps;
  n=0;
  while(n < steps)
    {
    char buf[100];
    double phi2=n*dphi;
    phi2-=PI;
    double x=phi2/(2.0*PI)*360;
    if(n==0)
      s.IMphi[n]=s.Mphi[n];
    else
      s.IMphi[n]=s.IMphi[n-1]+s.Mphi[n];
    sprintf(buf,"%g %g",x, (float)s.IMphi[n]);
    plot(buf);
    ++n;
    }
  plot( "e");

  fflush(plotpipe);
  plot("exit\n");
  fflush(plotpipe);
  fclose(plotpipe);
  }


 
Das ist meine persönliche Meinung dazu. Basierend auf einer nach bestem Wissen und Gewissen recherchierten Faktenlage.
asterix
Beiträge: 241
Registriert: 23.05.2018, 08:24

30.08.2019, 15:03

Die Ausgabe von Plots mit mehreren Y-Achsen macht in den üblichen Tabellenkalkulationsprogrammen immer Probleme, wenn man mehr als zwei Y-Achsen braucht. Das geht dort in der Regel nur, wenn man mehrere Diagramme trickreich übereinander legt. In Gnuplot funktioniert das deutlich einfacher. Man kann da beliebig viele Y-Achsen kreieren. Ich hab das einmal mit den über Open Office schon einmal dargestellten Klimadaten ausprobiert.

Das Ergebnis sieht dann so aus:

Abb1.gif
Abb1.gif (32.68 KiB) 122 mal betrachtet

Quelle_u_Daten.zip
(23.04 KiB) 3-mal heruntergeladen
Das ist meine persönliche Meinung dazu. Basierend auf einer nach bestem Wissen und Gewissen recherchierten Faktenlage.
asterix
Beiträge: 241
Registriert: 23.05.2018, 08:24

01.09.2019, 15:21

Jetzt noch die 3D-Darstellung. Da geht auch fast alles zunächst einmal von selbst. In Anlehnung zu der laufenden Klimadiskussion habe ich mir zum Ausprobieren einmal den Lorenz-Attraktor ausgesucht.  Ein Attraktor repräsentiert einen Phasenraum zu einer Problemlösung. Hier ein Strömungsproblem, das auch beim Wetter und beim Klima eine wesentliche Rolle spielt. Mit Phase ist dabei ein Zustand gemeint. Edward N. Lorenz hat zur Idealisierung eines hydrodynamischen Systems 3 Gleichungen entwickelt, über die sich alle möglichen Systemzustände zu einer Anordnung mit einem Fluid iterieren lassen. Das von ihm entwickelte System verhält sich chaotisch. Bei den gewählten Parametern deterministisch chaotisch. Die Realisierung ist ziemlich einfach. Um sich ein Verständnis zum Verhalten von chaotischen Systemen anzueignen, ist der Lorenz-Attraktor zwar nicht der geeignetste, weil er auf einer ziemlich komplexen Strömungslehre basiert, man kann aber ein paar grundsätzliche Sachen gut damit aufzeigen. Es gibt auch wesentlich allgemeiner gehaltene  Attraktoren. Ich habe gesehen, dass inzwischen sogar in einigen Schulen ein Verständnis zu chaotischen Systemen vermittelt wird. Einige grundsätzlichen Eigenschaften von Chaos lassen sich mit dem Lorenz-Attraktor aber auch recht einfach und ziemlich dekorativ aufzeigen, ohne dass man etwas von Strömungslehre verstehen muss.


Die Lorenz-Gleichungen beschreiben, wie sich ein Systemzustand nach einem winzigen Zeitabschnitt ändert. Die Parameter a.b und c beschreiben die Eigenschaften eines Fluids in einem System. X,y und z beschreiben dessen Systemzustand zu einer bestimmten Zeit.


    dx/ dt =  a * (y0 - x0)
    dy/ dt = (x0 * (b - z0) - y0)
    dz/ dt = (x0 * y0 - c * z0)


Der neue Systemzustand ergibt sich dann einfach aus den Addition von altem Zustand mit der errechneten Änderung.

 
Abb_1.gif
Abb_1.gif (8.86 KiB) 110 mal betrachtet

Das Programm ist ziemlich simpel:

Code: Alles auswählen

#include "stdio.h"
#include "stdlib.h"

FILE *plotpipe;
int open();
void plot(const char *x)
  {
  fprintf(plotpipe,"%s\n",x);
  }

main()
  {
  char buf[300];
  double xs=0,ys=0,zs=0;
  if(!(plotpipe = popen("gnuplot -persist","w")))
    {
    printf("Error opening pipe to GNU plot. Check if you have it! \n");
    exit(0);
    }
   int i=0;
   int steps=500000;
   double x0,y0,z0,x1,y1,z1;
   double dt = 0.0001;
   double a = 10.0;
   double b = 28.0;
   double c = 8.0/3.0;

   int mark_zustand=steps/2; // steps nach dem ein Zustand xs,ys,zs farblich markiert wird;
   int len_mark=1000;        // "länge" der Markierung

   x0 = 1.00;
   y0 = 1.0;
   z0 = 1.001;

   plot("set terminal wxt size 1200,600 enhanced font 'Verdana,10' persist");
   plot("set grid");
   plot("set key off"); // Legende aus
   plot("set xlabel \"x\"");
   plot("set ylabel \"y\"");
   plot("set zlabel \"z\"");
   plot("set ticslevel 0");
   plot("set style line 1 lt rgb \"blue\" lw 1");
   plot("set style line 2 lt rgb \"red\" lw 1");
   plot("set style line 3 lt rgb \"blue\" lw 4");
   plot("set style line 4 lt rgb \"red\" lw 1");
   plot("set style data  lines");
   sprintf(buf,"set label 1 \" %g,%g,%g\" at %g,%g,%g  nopoint  tc def",(float)x0,(float)y0,(float)z0,(float)x0,(float)y0,(float)z0);
   plot(buf);

   plot("splot '-'   using 1:2:3  ls 1, '-'   using 1:2:3  ls 2, '-'   using 1:2:3  ls 3, '-'   using 1:2:3  ls 4");       // ls 1 == style line 1
   sprintf(buf,"%g %g %g",(float)x0,(float)y0,(float)z0);
   plot(buf);
   for (i=0;i<steps;i++)
     {
     if(i == len_mark)               // weiter mit style line 2
       plot( "e");
     if(i == mark_zustand)           // weiter mit style line 3
       plot( "e");
     if(i == mark_zustand+len_mark)  // weiter mit style line 4
       {
       plot( "e");
       xs=x0;ys=y0;zs=z0;
       }
     x1 = x0 + dt * a * (y0 - x0);
     y1 = y0 + dt * (x0 * (b - z0) - y0);
     z1 = z0 + dt * (x0 * y0 - c * z0);
     x0 = x1;
     y0 = y1;
     z0 = z1;
     sprintf(buf,"%g %g %g",(float)x0,(float)y0,(float)z0);
     plot(buf);
     }
    plot( "e");
    if(mark_zustand > 0)
      {
      sprintf(buf,"set label 2 \"%.1f,%.1f,%.1f\"  at %g,%g,%g  nopoint  tc def",xs,ys,zs,(float)xs,(float)ys,(float)zs);
      plot(buf);
      plot( "replot");
      }
    fclose(plotpipe);
    return 0;
  }

Abb_2.gif
Abb_2.gif (90.76 KiB) 110 mal betrachtet
Das Ergebnis mit den Startwerten x0=1,y0=0,z0=0.

Jeder Punkt beschreibt einen Zustand. Die Kurven beschreiben sozusagen ein Gebiet mit Zuständen, in welches das durch die Parameter festgelegte System immer hinein will um dort dauerhaft gefangen zu bleiben.

Mit einem fernen Startpunkt kann man das verdeutlichen.

Abb_3.gif
Abb_3.gif (41.93 KiB) 110 mal betrachtet
 
Das Ergebnis mit den Startwerten x0=1,y0=-100,z0=0.

Trotz Berechenbarkeit ist aber keine Langzeitprognose möglich. Lorenz fand heraus, dass kleinste Änderungen bei den Startwerten nach einer festen Zeit zu völlig unterschiedlichen Sytemzuständen führen können. Das kann man auch leicht ausprobieren.

Abb_5.gif
Abb_5.gif (94.82 KiB) 110 mal betrachtet
 
Ergebnisse mit x0=1,y0=1,z0=1.001 und x0=1,y0=1,z0=1.002.
Der blaue Abschnitt markiert den Systemzustand nach 250000 Iterationen bei fast gleichem Startzustand. Man sieht den Grafiken schon an, dass die Bahnen trotz fast gleicher Startwerte ziemlich unterschiedlich verlaufen.

Die Sensitivität ist ein wichtiges Merkmal von Chaos. Kleinste Ungenauigkeiten verursachen mit der Zeit große Wirkungen. 250000 Iterationen entsprechen bei dt=0.0001s 25 Sekunden. Die 25 Sekunden erreicht man auch nach 125000 Iterationen bei dt=0.0002s. Die Genauigkeit bei der Integration geht dabei allerdings zurück.

Abb_6.gif
Abb_6.gif (61.15 KiB) 110 mal betrachtet
Das Ergebnis dazu wieder mit x0=1,y0=1,z0=1.001. Da kommt nach den vorgegebenen 25 Sekunden wieder was anderes raus.

Der Knackpunkt ist aber, dass der Bereich in dem die Systemzustände gefangen sind, trotz aller Widrigkeiten immer gleich bleibt. Aber nur so lange man die Parameter, welche die Eigenschaften des Systems festlegen, nicht ändert.



 
Das ist meine persönliche Meinung dazu. Basierend auf einer nach bestem Wissen und Gewissen recherchierten Faktenlage.
Antworten
  • Information