Finde kein Konvergenzkriterium für eine Iteration

[Bjoerk]

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;

[BUG]

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.

[Ultimator]

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.

[frankyboy1974]

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

[Bjoerk]

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.

[Bjoerk]

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;

[Mikkey]

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.

[Bjoerk]

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;