線形方程式を解く

Question

C、Objective C、または（必要な場合）C++の線形方程式のシステムをプログラムで解決する必要があります。私はa、b、およびtxための最良の近似値を取得したいのですが、このことから

-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 

：

は、ここでの式の例です。

Answer 1

Cramer's Rule と Gaussian Eliminationは、二つの良い、汎用アルゴリズム（またSimultaneous Linear Equationsを参照）です。コードをお探しの場合は、GiNaC、Maxima、SymbolicC++（ライセンス要件に応じて異なります）をご覧ください。

編集：私はあなたがC言語で作業していることを知っていますが、SymPy（Pythonでコンピュータ代数システム）の良い言葉を入れなければなりません。あなたはそのアルゴリズムから多くのことを学ぶことができます（あなたが少しのPythonを読むことができるなら）。また、新しいBSDライセンスの下にありますが、無料の数学パッケージのほとんどはGPLです。

Answer 2

作業を行うか、実際にマトリックス操作などを行い、各ステップを実行するソフトウェアパッケージをお探しですか？

最初に、私の同僚がちょうどOcaml GLPKを使用しました。それはちょうどGLPKのための包みですが、それは物事を設定する多くのステップを削除します。あなたはGLPKに固執しなければならないようですが、Cでもそうです。後者については、私は以前にLPを学んでいた古い記事を保存するためのおいしいおかげで、PDF。具体的なヘルプ設定が必要な場合は、私たちに知らせてください。私や誰かがここに戻り、助けてくれるでしょうが、私はそれはかなり正直だと思います。がんばろう！

Answer 3

1次方程式の3x3システムでは、独自のアルゴリズムを公開することは大丈夫でしょう。

ただし、精度、ゼロ除算、実際に小さな数値、無限に多くのソリューションについて何をするかについて心配する必要があります。私の提案は、のような標準数値線形代数パッケージを使うことです。

Answer 4

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

この本は、なぜアルゴリズムがうまくいくのかを教え、それらのアルゴリズムのかなりうまくデバッグされた実装を示します。

もちろん、盲目的にはCLAPACKを使用しましたが（私はこれを大成功に使用しました）、最初にGaussian Eliminationアルゴリズムを手入力して、少なくとも行った作業これらのアルゴリズムを安定させることになる。

もっと興味深い線形代数を使っているなら、Octaveのソースコードを見てみると、たくさんの質問に答えることができます。

Answer 5

Template Numerical Toolkit NISTのNISTには、このようなツールがあります。

信頼性の高い方法の1つは、QR Decompositionを使用することです。

ここでは、私のコードで "GetInverse（A、InvA）"を呼び出すことができるようにラッパーの例を示し、インバースにインバースを入れます。

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

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

Answer 6

あなたの質問の文言から、あなたは未知数より多くの方程式を持っているようで、不一致を最小限に抑えたいと思っています。これは、典型的には、線形回帰を用いて行われ、不一致の二乗の和を最小化する。データのサイズに応じて、スプレッドシートまたは統計パッケージで行うことができます。 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 

Answer 7

、あなたは常に変数と方程式の数が同じになる場合は他の人がI.より良い答えていることが、実装が簡単だとして、私はCramer's ruleが好き。行列の行列式を計算する関数を書いてください（またはすでに書かれているものを使用してください）。そして、2つの行列の行列式を分割します。

Answer 8

プログラムでこれを解決する方法は、手で（掛け算と減算を行い、結果を方程式に戻す）まったく同じ方法で解決できます。これはかなり標準的な中等学校レベルの数学です。

-44.3940 = 50a + 37b + c (A) 
-45.3049 = 43a + 39b + c (B) 
-44.9594 = 52a + 41b + c (C) 

(A-B): 0.9109 = 7a - 2b (D) 
(B-C): 0.3455 = -9a - 2b (E) 

(D-E): 1.2564 = 16a (F) 

(F/16): a = 0.078525 (G) 

Feed G into D: 
     0.9109 = 7a - 2b 
    => 0.9109 = 0.549675 - 2b (substitute a) 
    => 0.361225 = -2b (subtract 0.549675 from both sides) 
    => -0.1806125 = b (divide both sides by -2) (H) 

Feed H/G into A: 
     -44.3940 = 50a + 37b + c 
    => -44.3940 = 3.92625 - 6.6826625 + c (substitute a/b) 
    => -41.6375875 = c (subtract 3.92625 - 6.6826625 from both sides) 

だからで終わる：あなたが戻ってA、BおよびCにこれらの値を接続した場合、あなたは彼らが正しいです見つける

a = 0.0785250 
b = -0.1806125 
c = -41.6375875 

。

トリックは、単純に3x2マトリックスに減らした単純な4x3マトリックスを使用し、次に2x1を "a = n"にします.nは実際の数です。いったんそれを持っていれば、それを次の行列に入れて別の値にし、次にそれらの2つの値をすべての変数を解くまで次の行列に入れます。

N個の異なる方程式があれば、N個の変数を常に解くことができます。これら二つではないので、私は明確な言う：

7a + 2b = 50 
14a + 4b = 100 

彼らは2倍同じ式ですので、あなたは彼らから解決策を得ることができない - あなたは真のが、役に立たない文で2枚のその後、減算の葉で最初に掛けます：一例として

0 = 0 + 0 

---

、ここでは、あなたの質問に置かれている連立方程式をうまくいく、いくつかのCコードです。最初のいくつかの必要なタイプ、変数、数式を印刷するための支持機能、およびmain開始：

#include <stdio.h> 

typedef struct { double r, a, b, c; } tEquation; 
tEquation equ1[] = { 
    { -44.3940, 50, 37, 1 },  // -44.3940 = 50a + 37b + c (A) 
    { -45.3049, 43, 39, 1 },  // -45.3049 = 43a + 39b + c (B) 
    { -44.9594, 52, 41, 1 },  // -44.9594 = 52a + 41b + c (C) 
}; 
tEquation equ2[2], equ3[1]; 

static void dumpEqu (char *desc, tEquation *e, char *post) { 
    printf ("%10s: %12.8lf = %12.8lfa + %12.8lfb + %12.8lfc (%s)\n", 
     desc, e->r, e->a, e->b, e->c, post); 
} 

int main (void) { 
    double a, b, c; 

次に二つの未知数を有する2つの方程式の3つの未知数を有する3つの式の還元：

// First step, populate equ2 based on removing c from equ. 

    dumpEqu (">", &(equ1[0]), "A"); 
    dumpEqu (">", &(equ1[1]), "B"); 
    dumpEqu (">", &(equ1[2]), "C"); 
    puts (""); 

    // A - B 
    equ2[0].r = equ1[0].r * equ1[1].c - equ1[1].r * equ1[0].c; 
    equ2[0].a = equ1[0].a * equ1[1].c - equ1[1].a * equ1[0].c; 
    equ2[0].b = equ1[0].b * equ1[1].c - equ1[1].b * equ1[0].c; 
    equ2[0].c = 0; 

    // B - C 
    equ2[1].r = equ1[1].r * equ1[2].c - equ1[2].r * equ1[1].c; 
    equ2[1].a = equ1[1].a * equ1[2].c - equ1[2].a * equ1[1].c; 
    equ2[1].b = equ1[1].b * equ1[2].c - equ1[2].b * equ1[1].c; 
    equ2[1].c = 0; 

    dumpEqu ("A-B", &(equ2[0]), "D"); 
    dumpEqu ("B-C", &(equ2[1]), "E"); 
    puts (""); 

次に、未知の1と1つの方程式には2つの未知数を有する2次方程式の削減：我々はの式を持っていることを今

// Next step, populate equ3 based on removing b from equ2. 

    // D - E 
    equ3[0].r = equ2[0].r * equ2[1].b - equ2[1].r * equ2[0].b; 
    equ3[0].a = equ2[0].a * equ2[1].b - equ2[1].a * equ2[0].b; 
    equ3[0].b = 0; 
    equ3[0].c = 0; 

    dumpEqu ("D-E", &(equ3[0]), "F"); 
    puts (""); 

タイプnumber1 = unknown * number2の場合、unknown <- number1/number2で未知の値を計算することができます。次に、その値を計算したら、2つの未知数を持つ方程式の1つに代入し、2番目の値を計算します。その後、元の方程式の一つに両方のものを（現在知られている）の未知数を代用し、あなたが今、すべての3つの未知数の値を持っている：

：

// Finally, substitute values back into equations. 

    a = equ3[0].r/equ3[0].a; 
    printf ("From (F ), a = %12.8lf (G)\n", a); 

    b = (equ2[0].r - equ2[0].a * a)/equ2[0].b; 
    printf ("From (D,G ), b = %12.8lf (H)\n", b); 

    c = (equ1[0].r - equ1[0].a * a - equ1[0].b * b)/equ1[0].c; 
    printf ("From (A,G,H), c = %12.8lf (I)\n", c); 

    return 0; 
} 

そのコードの出力は、この答えでは、以前の計算と一致しました

  >: -44.39400000 = 50.00000000a + 37.00000000b + 1.00000000c (A) 
     >: -45.30490000 = 43.00000000a + 39.00000000b + 1.00000000c (B) 
     >: -44.95940000 = 52.00000000a + 41.00000000b + 1.00000000c (C) 

     A-B: 0.91090000 = 7.00000000a + -2.00000000b + 0.00000000c (D) 
     B-C: -0.34550000 = -9.00000000a + -2.00000000b + 0.00000000c (E) 

     D-E: -2.51280000 = -32.00000000a + 0.00000000b + 0.00000000c (F) 

From (F ), a = 0.07852500 (G) 
From (D,G ), b = -0.18061250 (H) 
From (A,G,H), c = -41.63758750 (I) 

Answer 9

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

が出力されます。
===ソルバー財団サービスレポート===
日時：2009年4月20日23： 29:55
モデル名：要求されたデフォルト
機能：LP
解決時間（ミリ秒）：1027
合計時間（ミリ秒）：1414
最適
ソルバー選択：Microsoft.SolverFoundation.Solvers.SimplexSolver
ディレクティブ：
Microsoft.SolverFoundation.Services.Directive
アルゴリズム：プライマル
算術：ハイブリッド
価格（正確に）：デフォルト完了ステータスを解決 価格（ダブル）：SteepestEdge
根拠：3
===ソリューションの詳細===
目標：
カウント
ピボットスラックの 決定：
：0.0785250000000004
B：-0.180612500000001
C：-41.6375875

Answer 10

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