Finde kein Konvergenzkriterium für eine Iteration

[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;