AGB  ·  Datenschutz  ·  Impressum  







Anmelden
Nützliche Links
Registrieren
Zurück Delphi-PRAXiS Sprachen und Entwicklungsumgebungen Object-Pascal / Delphi-Language Finde kein Konvergenzkriterium für eine Iteration

Finde kein Konvergenzkriterium für eine Iteration

Ein Thema von Bjoerk · begonnen am 10. Nov 2015 · letzter Beitrag vom 12. Nov 2015
Antwort Antwort
Bjoerk

Registriert seit: 28. Feb 2011
Ort: Mannheim
1.376 Beiträge
 
Delphi 2007 Professional
 
#1

Finde kein Konvergenzkriterium für eine Iteration

  Alt 10. Nov 2015, 20:39
Delphi-Version: 2007
Ich habe eine Iteration. Ausgehend von einem Wert kX wir ein erster Wert "Force" ermittelt. Anschließend werden 3 weitere Werte ermittelt. Das geschieht aus 3 Gleichungen mit 3 Unbekannten (-> "TopRight", "BottomLeft" und "BottomRight"). Diese Werte werden in die Gleichungen eingesetzt und überprüft, ob alle 3 Gleichungen Null werden. Ist das so, dann hat man das richtige kX gefundent. Mein Problem hier ist, daß ich kein Konvergenzkriterium finde. Es iteriert nicht. Die Werte springen. Wie könnte man denn das machen? Probiere schon ein ganze Weile..

Delphi-Quellcode:
function TBoltApproximationEx.CalcDefault: boolean; // KX = 0 .. 1;
var
  I: integer;
  OldValue: double;
  NewValue: double;
  Alpha: double;
  KxV: double;
  dKxV: double;
  XV: double;
  C: double;
  B: double;
  D: double;
  X: double;
  hx: double;
  hy: double;
  c1: double;
  d1: double;
  aX: double;
  aY: double;
  TopRight: double;
  BottomLeft: double;
  BottomRight: double;
  Force: double;
  SigmaD: double;
  Mx: double;
  My: double;
  SumN: double;
  SumMx: double;
  SumMy: double;
  TopRightX: double;
  TopRightY: double;
  BottomLeftX: double;
  BottomLeftY: double;
begin
  FStrings.Clear;

  Alpha := GetAlpha;
  C := Cos(Alpha);

  Mx := Abs(FMX);
  My := Abs(FMY);

  b := FB / 1000; // m;
  d := FD / 1000;
  c1 := FC1 / 1000;
  d1 := FD1 / 1000;

  hx := d - c1; // c2 = c1;
  hy := b - d1; // d2 = d1;

  TopRightX := b - d1; // TopRight;
  TopRightY := -c1;
  BottomLeftX := d1; // BottomLeft;
  BottomLeftY := -d + c1;

  KxV := 0;
  dKxV := 0.001;

  NewValue := 0;
  I := 0;
  repeat // KxV;
    Inc(I);
    OldValue := NewValue;

    XV := KxV * hx; // m;
    FPolygon.SetFromBemeRect(b, d, XV, Alpha);
    FPolygon.Calc;

    X := XV * C;
    SigmaD := 0.5 * Min(FFy * 1E3, FFy * 1E3 * X / (hx - X)); // kN/m2; FFy = Material N/mm2;

    Force := -FPolygon.Area * SigmaD; // kN;
    aX := FPolygon.CenterX; // m;
    aY := -FPolygon.CenterY; // m;

    if FPolygon.PtIn(TopRightX, TopRightY) and FPolygon.PtIn(BottomLeftX, BottomLeftY) then
    begin
      TopRight := 0;
      BottomLeft := 0;
      BottomRight := FN - Force;
    end
    else
      if FPolygon.PtIn(TopRightX, TopRightY) then
      begin
        TopRight := 0;

        FGauss.Count := 2;

        FGauss.A[0, 0] := hx;
        FGauss.A[0, 1] := hx;
        FGauss.B[0] := Mx + FN * d / 2 - Force * aY;

        FGauss.A[1, 0] := d1;
        FGauss.A[1, 1] := hy;
        FGauss.B[1] := My + FN * b / 2 - Force * aX;

        FGauss.Complete;
        FGauss.Solve;
        BottomLeft := Max(FGauss.X[0], 0);
        BottomRight := Max(FGauss.X[1], 0);
      end
      else
        if FPolygon.PtIn(BottomLeftX, BottomLeftY) then
        begin
          BottomLeft := 0;

          FGauss.Count := 2;

          FGauss.A[0, 0] := hx;
          FGauss.A[0, 1] := c1;
          FGauss.B[0] := Mx + FN * d / 2 - Force * aY;

          FGauss.A[1, 0] := hy;
          FGauss.A[1, 1] := hy;
          FGauss.B[1] := My + FN * b / 2 - Force * aX;

          FGauss.Complete;
          FGauss.Solve;
          BottomRight := Max(FGauss.X[0], 0);
          TopRight := Max(FGauss.X[1], 0);
        end
        else
        begin
          FGauss.Count := 3;

          FGauss.A[0, 0] := hx;
          FGauss.A[0, 1] := hx;
          FGauss.A[0, 2] := c1;
          FGauss.B[0] := Mx + FN * d / 2 - Force * aY;

          FGauss.A[1, 0] := d1;
          FGauss.A[1, 1] := hy;
          FGauss.A[1, 2] := hy;
          FGauss.B[1] := My + FN * b / 2 - Force * aX;

          FGauss.A[2, 0] := 1;
          FGauss.A[2, 1] := 1;
          FGauss.A[2, 2] := 1;
          FGauss.B[2] := FN - Force;

          FGauss.Complete;
          FGauss.Solve;
          BottomLeft := Max(FGauss.X[0], 0);
          BottomRight := Max(FGauss.X[1], 0);
          TopRight := Max(FGauss.X[2], 0);
        end;

    // Bedingung 1: Summe Mx = 0;
    SumMx := -Mx -FN * d / 2 + Force * aY + BottomLeft * hx + BottomRight * hx + TopRight * c1; // kNm;
    // Bedingung 2: Summe My = 0;
    SumMy := -My -FN * b / 2 + Force * aX + BottomLeft * d1 + BottomRight * hy + TopRight * hy; // kNm;
    // Bedingung 3: Summe N = 0;
    SumN := -FN + Force + BottomLeft + BottomRight + TopRight; // kN;

    NewValue := SumN; // oder SumMx oder SumMy ???
    if NewValue * OldValue < 0 then dKxV := -dKxV / 2;

    KxV := KxV + dKxV;

    FStrings.Add(Format('I %4d N %4d X %.2f SumMx %.4f SumMy %.4f SumN %.4f Force %.4f TopRight %.4f BottomLeft %.4f BottomRight %.4f',
      [I, FGauss.Count, 100 * X, SumMx, SumMy, SumN, Force, TopRight, BottomLeft, BottomRight]));

  until (I = 10000) or (Abs(NewValue) < 0.01);
end;
  Mit Zitat antworten Zitat
Benutzerbild von BUG
BUG

Registriert seit: 4. Dez 2003
Ort: Cottbus
2.094 Beiträge
 
#2

AW: Finde kein Konvergenzkriterium für eine Iteration

  Alt 10. Nov 2015, 21:46
Die traurige Wahrheit ist: nicht alle Iterationen konvergieren.
Wenn du etwas mehr erzählst was du da eigendlich machst, fällt vielleicht jemandem was ein.
  Mit Zitat antworten Zitat
Benutzerbild von Ultimator
Ultimator

Registriert seit: 17. Feb 2004
Ort: Coburg
1.860 Beiträge
 
FreePascal / Lazarus
 
#3

AW: Finde kein Konvergenzkriterium für eine Iteration

  Alt 11. Nov 2015, 07:11
Ich hab mir deinen Code jetzt nicht im Detail angeschaut, aber so Sachen wie N/mm² und kX klingen für mich nach einem Finite-Elemente-/Finite-Differenzen-Code.

Bei solchen iterativen Geschichten kommt es immer auf die Schrittweite an. Ist die zu groß, ist der Code instabil. Vielleicht wäre es hilfreich, wenn du deinen berechneten Wert für Verschiebungen und Kräfte nicht einfach für die nächste Iteration übernimmst, sondern mit einem Faktor "unterrelaxierst":

Wert_fuer_naechste_Iteration := alpha * Wert_der_aktuellen_Iteration + (1 - alpha) * Soeben_ausgerechneter_Wert Mit alpha kannst du dann so lange spielen, bis du einen guten Kompromiss zwischen Stabilität und Geschwindigkeit gefunden hast.
Julian J. Pracht
  Mit Zitat antworten Zitat
Benutzerbild von frankyboy1974
frankyboy1974

Registriert seit: 7. Apr 2015
Ort: SH
169 Beiträge
 
Delphi XE7 Professional
 
#4

AW: Finde kein Konvergenzkriterium für eine Iteration

  Alt 11. Nov 2015, 07:46
hallo,

ohne das eigentliche Problem zu verstehen, würde ich das Problem in der Verwendung der Min/Max-Funktionen vermuten, daraus ergibt sich dann eine unstetige Funktion, die zwischen zwei Werten hin-und her springt und eben nicht konvergiert.

mfg
Java ist auch eine Insel.
Ist Delphi von Oracle?
In meiner Buchstabensuppen fehlt das C++!
  Mit Zitat antworten Zitat
Bjoerk

Registriert seit: 28. Feb 2011
Ort: Mannheim
1.376 Beiträge
 
Delphi 2007 Professional
 
#5

AW: Finde kein Konvergenzkriterium für eine Iteration

  Alt 11. Nov 2015, 08:45
Danke daß ihr euch für mein Problem interessiert. Wenn ich da aber anfange zu erläutern, wird‘s ein Baustatik Einsteigerseminar? Ich zieh das ganze mal etwas anders auf. Ich denke, man kann Unbekannte zu Resultierenden zusammenfassen. Vielleicht iteriert es dann besser. Mal sehen.. Melde mich ggf. nochmal. Dann kann auch das Max entfallen. Für ein ebenes Problem ist es hier ganz nett erläutert.
Angehängte Dateien
Dateityp: pdf Kopfplattenanschluss.pdf (101,7 KB, 27x aufgerufen)
  Mit Zitat antworten Zitat
Bjoerk

Registriert seit: 28. Feb 2011
Ort: Mannheim
1.376 Beiträge
 
Delphi 2007 Professional
 
#6

AW: Finde kein Konvergenzkriterium für eine Iteration

  Alt 11. Nov 2015, 20:06
Jo, so wie ich das vorhatte, jeeht das net. Habs mir mal aufskizziert. Ist anders als im Stahlbetonbau, die Kräfte ergeben sich hier rein aus der Spannungsverteilung. 08/15 Statik müßte man halt können ..

Code:
Z.i (i = 0..3) = Abs(Sigma.D) * (h.i – x.v) / x.v / Sum.Sigma.Z * Abs(Force), für alle h.i > x.v
Code:
repeat
  Alpha := ..
  repeat
    XV := ..
  until SumMx = 0;
until SumMy = 0;
Angehängte Dateien
Dateityp: pdf KopfplattenanschlussNeu.pdf (1,51 MB, 21x aufgerufen)
  Mit Zitat antworten Zitat
Mikkey

Registriert seit: 5. Aug 2013
265 Beiträge
 
#7

AW: Finde kein Konvergenzkriterium für eine Iteration

  Alt 12. Nov 2015, 10:06
Ohne Dir jetzt nahetreten zu wollen, aber die Zeichnung macht ziemlich genau denselben Eindruck wie der Code in Deiner Frage.

Was hältst Du davon, beides ein wenig zu ordnen (um Deiner selbst willen)?

Gib die Aufgabenstellung an, ich bin zwar kein Statiker, aber ich habe innerhalb der angewandten Mathematik auch solche Aufgaben lösen müssen. Deshalb solltest Du die Herleitung des Algorithmus aus der Differentialgleichung und den Rand/Anfangswerten angeben.

Teile die ellenlange Funktion mit den gefühlten 200 Variablen in Funktionen auf, die einzelne Aufgaben erledigen können. Dann werden sich auch mehr Teilnehmer hier den Code zu Gemüte führen.
  Mit Zitat antworten Zitat
Bjoerk

Registriert seit: 28. Feb 2011
Ort: Mannheim
1.376 Beiträge
 
Delphi 2007 Professional
 
#8

AW: Finde kein Konvergenzkriterium für eine Iteration

  Alt 12. Nov 2015, 13:50
Glaub ich dir, glaub ich dir. Und ja, der Code hatte noch Potential (Siehe unten). Die Skizze sollte auch keinen Schönheitswettbewerb gewinnen. Was ist dir noch unklar? Vielleicht nochmal zum Procedere: Wir raten ein Alpha und wir raten ein xV. Ob wir richtig geraten haben, können wir überprüfen. Auf der Höhe xV schneidet eine Senkrechte dieser Line den Querschnitt b/d. Wir erhalten eine Fläche. Diese Fläche ist die Grundfläche eines Polyeders. An den Schnittpunkten ist Sigma Null. In der linken oberen Ecke maximal. Die anderen Eckpunkte kann man linear interpolieren. Das Volumen dieses Polyeders steht stellvertretend für die resultierende Druckkraft im Querschnitt. Diese brauchen wir um zu überprüfen ob wir richtig geraten haben.

Delphi-Quellcode:
function TBoltApproximationD3.Calc: boolean; // KX = 0 .. 1;
const
  cnst = 1000;
  cnstdKxV = 0.001;
  cnstdAlpha = Pi / 180;
  cEps = 0.0001;
  cTopLeft = 0;
  cTopRight = 3;
  cBottomLeft = 1;
  cBottomRight = 2;
var
  Exchanged: boolean;
  Alpha, dAlpha, KxV, dKxV, xV, SumMx, SumMy, ForceCenterX, Temp: double;
begin
  Result := false;
  ResultClear;
  if CanCalc then
  begin
    Exchanged := false;
    if GetAlpha > Pi / 4 then
    begin
      Exchanged := true;
      Exchange(FMx, FMy);
      Exchange(FB, FD);
      Exchange(FC1, FD1);
    end;
    try
      FPoints.Fyd := FFy;
      // Alpha;
      dAlpha := cnstdAlpha;
      Alpha := -dAlpha;
      SumMy := 0;
      repeat
        Alpha := Alpha + dAlpha;
        FPoints.Alpha := Alpha; // Rad;
        FPoints.Count := 4;
        FPoints.X[cTopLeft] := FD1 / cnst;
        FPoints.Y[cTopLeft] := -FC1 / cnst;
        FPoints.X[cBottomLeft] := FD1 / cnst;
        FPoints.Y[cBottomLeft] := -FD / cnst + FC1 / cnst;
        FPoints.X[cBottomRight] := FB / cnst - FD1 / cnst;
        FPoints.Y[cBottomRight] := -FD / cnst + FC1 / cnst;
        FPoints.X[cTopRight] := FB / cnst - FD1 / cnst;
        FPoints.Y[cTopRight] := -FC1 / cnst;
        // KxV;
        dKxV := cnstdKxV;
        KxV := -dKxV;
        SumMx := 0;
        repeat
          KxV := KxV + dKxV;
          xV := KxV * FPoints.MaxH; // m;
          FPoints.xV := xV;
          if FUsePolyeder then
          begin
            FPoints.SigmaD := Min(FFy * cnst, FFy * cnst * xV / (FPoints.MaxH - xV));
            FPoints.Force := FPolyeder.SolidStressFromRect(FB / cnst, FD / cnst, xV, FPoints.SigmaD, Alpha);
            ForceCenterX := FPolyeder.Center.X;
          end
          else
          begin
            FPoints.SigmaD := 0.5 * Min(FFy * cnst, FFy * cnst * xV / (FPoints.MaxH - xV));
            FPoints.Force := FPolygon.SimplifiedStressFromRect(FB / cnst, FD / cnst, xV, FPoints.SigmaD, Alpha);
            ForceCenterX := FPolygon.CenterX;
          end;
          Temp := Abs(FMx) - FPoints.SumMx;
          if SumMx * Temp < 0 then
            dKxV := -dKxV / 2;
          SumMx := Temp;
          Result := (CompareValue(KxV, 0) >= 0) and (CompareValue(KxV, 1) <= 0);
        until (Abs(SumMx) < cEps) or (not Result);
        if Result then
        begin
          Temp := Abs(FMy) - FPoints.SumMy + FPoints.Force * ForceCenterX;
          if SumMy * Temp < 0 then
            dAlpha := -dAlpha / 2;
          SumMy := Temp;
        end;
        Result := (CompareValue(Alpha, 0) >= 0) and (CompareValue(Alpha, Pi / 2) <= 0);
      until (Abs(SumMy) < cEps) or (not Result);
    finally
      if Exchanged then
      begin
        Exchange(FMx, FMy);
        Exchange(FB, FD);
        Exchange(FC1, FD1);
      end;
      if not Result then
        ResultClear;
    end;
  end;
end;

Geändert von Bjoerk (12. Nov 2015 um 13:55 Uhr)
  Mit Zitat antworten Zitat
Themen-Optionen Thema durchsuchen
Thema durchsuchen:

Erweiterte Suche
Ansicht

Forumregeln

Es ist dir nicht erlaubt, neue Themen zu verfassen.
Es ist dir nicht erlaubt, auf Beiträge zu antworten.
Es ist dir nicht erlaubt, Anhänge hochzuladen.
Es ist dir nicht erlaubt, deine Beiträge zu bearbeiten.

BB-Code ist an.
Smileys sind an.
[IMG] Code ist an.
HTML-Code ist aus.
Trackbacks are an
Pingbacks are an
Refbacks are aus

Gehe zu:

Impressum · AGB · Datenschutz · Nach oben
Alle Zeitangaben in WEZ +1. Es ist jetzt 10:26 Uhr.
Powered by vBulletin® Copyright ©2000 - 2019, Jelsoft Enterprises Ltd.
LinkBacks Enabled by vBSEO © 2011, Crawlability, Inc.
Delphi-PRAXiS (c) 2002 - 2019 by Daniel R. Wolf