48

C、Objective C、または（必要に応じて）C++で連立一次方程式をプログラムで解く必要があります。

方程式の例を次に示します。

-44.3940 = a * 50.0 + b * 37.0 + tx
-45.3049 = a * 43.0 + b * 39.0 + tx
-44.9594 = a * 52.0 + b * 41.0 + tx

aこのことから、、、、の最適な近似値を取得したいと思いbますtx。

4

11 に答える 11

19

クラメルの法則 と ガウスの消去法 は、2つの優れた汎用アルゴリズムです（連立一次方程式も参照）。コードをお探しの場合は、GiNaC、Maxima、SymbolicC ++をチェックしてください（もちろん、ライセンス要件によって異なります）。

編集：あなたがCランドで働いていることは知っていますが、SymPy（Pythonの数式処理システム）についても良い言葉を述べなければなりません。そのアルゴリズムから多くのことを学ぶことができます（Pythonを少し読むことができる場合）。また、それは新しいBSDライセンスの下にありますが、無料の数学パッケージのほとんどはGPLです。

于 2008-08-03T18:37:24.003 に答える
7

3x3 の線形方程式系の場合、独自のアルゴリズムを展開しても問題ないと思います。

ただし、精度、ゼロまたは非常に小さな数による除算、および無限に多くの解について何をすべきかについて心配する必要がある場合があります。私の提案は、 LAPACKなどの標準の数値線形代数パッケージを使用することです。

于 2008-08-08T17:59:27.907 に答える
7

Microsoft Solver Foundationをご覧ください。

これを使用すると、次のようなコードを記述できます。

  SolverContext context = SolverContext.GetContext();
  Model model = context.CreateModel();

  Decision a = new Decision(Domain.Real, "a");
  Decision b = new Decision(Domain.Real, "b");
  Decision c = new Decision(Domain.Real, "c");
  model.AddDecisions(a,b,c);
  model.AddConstraint("eqA", -44.3940 == 50*a + 37*b + c);
  model.AddConstraint("eqB", -45.3049 == 43*a + 39*b + c);
  model.AddConstraint("eqC", -44.9594 == 52*a + 41*b + c);
  Solution solution = context.Solve();
  string results = solution.GetReport().ToString();
  Console.WriteLine(results);

出力は次のとおりです:
===Solver Foundation Service Report===
Datetime: 04/20/2009 23:29:55
Model Name: Default
要求された機能: LP
Solve Time (ms): 1027
Total Time (ms): 1414
Solve完了ステータス: 最適な
ソルバーが選択されました: Microsoft.SolverFoundation.Solvers.SimplexSolver
ディレクティブ:
Microsoft.SolverFoundation.Services.Directive
アルゴリズム:主要な
算術演算: ハイブリッド
価格 (正確): デフォルト
の価格 (倍精度): SteepestEdge
基準: スラック
ピボット カウント: 3
== =ソリューションの詳細===
目標:

決定:
a: 0.0785250000000004
b: -0.180612500000001
c: -41.6375875

于 2009-04-21T03:33:42.160 に答える
3

作業を実行するソフトウェアパッケージを探していますか、それとも実際に行列演算などを実行して各ステップを実行しますか？

最初の、私の同僚はちょうどOcamlGLPKを使用しました。これはGLPKの単なるラッパーですが、設定の多くの手順が削除されます。ただし、CではGLPKを使用する必要があるようです。後者については、私がしばらく前にLPを学んでいた古い記事、 PDFを保存してくれたおいしいおかげで。さらにセットアップについて具体的なサポートが必要な場合は、お知らせください。私または誰かが戻ってきてサポートしてくれると確信していますが、ここからはかなり簡単だと思います。幸運を！

于 2008-08-03T18:27:13.180 に答える
3

NIST のTemplate Numerical Toolkitには、それを行うためのツールがあります。

より信頼できる方法の 1 つは、QR Decompositionを使用することです。

コードで「GetInverse(A, InvA)」を呼び出すことができるラッパーの例を次に示します。これにより、逆が InvA に配置されます。

void GetInverse(const Array2D<double>& A, Array2D<double>& invA)
   {
   QR<double> qr(A);  
   invA = qr.solve(I); 
   }

Array2D はライブラリで定義されています。

于 2008-08-25T18:53:32.837 に答える
3

実行時間の効率に関しては、他の人が私よりも良い回答をしています。常に変数と同じ数の方程式を使用する場合、実装が簡単なクラマーの規則が気に入っています。行列の行列式を計算する関数を作成するだけです (または、既に作成されている関数を使用してください。きっと見つかるはずです)。2 つの行列の行列式を除算します。

于 2008-09-19T03:36:13.207 に答える
2

個人的には、 Numerical Recipesのアルゴリズムが好きです。(私は C++ 版が好きです。)

この本では、アルゴリズムが機能する理由を説明し、さらに、これらのアルゴリズムのかなり適切にデバッグされた実装を示します。

もちろん、盲目的にCLAPACKを使用することもできます(私はこれを使用して大成功を収めました) が、少なくともこれらのアルゴリズムの作成に費やされた作業の種類についてかすかな考えを持つために、最初にガウス消去アルゴリズムを手動で入力します。安定。

後で、もっと興味深い線形代数をやっているなら、Octaveのソース コードを調べてみると、多くの疑問が解決します。

于 2008-08-25T18:22:48.010 に答える
2

質問の文言から、未知数よりも多くの方程式があり、矛盾を最小限に抑えたいと思われます。これは通常、不一致の二乗和を最小化する線形回帰で行われます。データのサイズに応じて、スプレッドシートまたは統計パッケージでこれを行うことができます。R は、線形回帰などを行う高品質の無料パッケージです。線形回帰には多くのこと (および多くの落とし穴) がありますが、単純なケースでは簡単に行うことができます。データを使用した R の例を次に示します。「tx」はモデルへのインターセプトであることに注意してください。

> y <- c(-44.394, -45.3049, -44.9594)
> a <- c(50.0, 43.0, 52.0)
> b <- c(37.0, 39.0, 41.0)
> regression = lm(y ~ a + b)
> regression

Call:
lm(formula = y ~ a + b)

Coefficients:
(Intercept)            a            b  
  -41.63759      0.07852     -0.18061
于 2008-09-17T14:22:56.880 に答える
1 
function x = LinSolve(A,y)
%
% Recursive Solution of Linear System Ax=y
% matlab equivalent: x = A\y 
% x = n x 1
% A = n x n
% y = n x 1
% Uses stack space extensively. Not efficient.
% C allows recursion, so convert it into C. 
% ----------------------------------------------
n=length(y);
x=zeros(n,1);
if(n>1)
    x(1:n-1,1) = LinSolve( A(1:n-1,1:n-1) - (A(1:n-1,n)*A(n,1:n-1))./A(n,n) , ...
                           y(1:n-1,1) - A(1:n-1,n).*(y(n,1)/A(n,n))); 
    x(n,1) = (y(n,1) - A(n,1:n-1)*x(1:n-1,1))./A(n,n); 
else
    x = y(1,1) / A(1,1);
end
于 2014-12-30T15:24:38.107 に答える