2008-08-03 96 views
34

मुझे सी, उद्देश्य सी, या (यदि आवश्यक हो) सी ++ में रैखिक समीकरणों की एक प्रणाली को प्रोग्रामेटिक रूप से हल करने की आवश्यकता है।एक रैखिक समीकरण को हल करना

यहाँ समीकरणों के एक उदाहरण है:

-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 के लिए सबसे अच्छा सन्निकटन प्राप्त करना चाहते हैं।

+0

अन्य लोगों को इस का जवाब दे दिया है, लेकिन पुस्तक * संख्यात्मक विश्लेषण की जाँच: वैज्ञानिक कम्प्यूटिंग के गणित * किनकैड द्वारा और चेनी। पुस्तक बड़े पैमाने पर समीकरणों की विभिन्न प्रणालियों को हल करने के बारे में है। – Matthew

उत्तर

17

Cramer's Rule और Gaussian Elimination दो अच्छा, सामान्य उद्देश्य एल्गोरिदम (भी Simultaneous Linear Equations देखें) कर रहे हैं। यदि आप कोड की तलाश में हैं, तो GiNaC, Maxima, और SymbolicC++ देखें (निश्चित रूप से आपकी लाइसेंसिंग आवश्यकताओं के आधार पर)।

संपादित करें: मुझे पता है कि आप सी भूमि में काम कर रहे हैं, लेकिन मुझे SymPy (पायथन में एक कंप्यूटर बीजगणित प्रणाली) के लिए एक अच्छा शब्द भी देना होगा। आप इसके एल्गोरिदम से बहुत कुछ सीख सकते हैं (यदि आप थोड़ा पाइथन पढ़ सकते हैं)। इसके अलावा, यह नए बीएसडी लाइसेंस के तहत है, जबकि अधिकांश मुफ्त गणित पैकेज जीपीएल हैं।

+12

असल में, न तो क्रैमर का शासन और न ही गाऊशियन उन्मूलन वास्तविक दुनिया में बहुत अच्छा है। न तो अच्छे संख्यात्मक गुण हैं, और न ही गंभीर अनुप्रयोगों के लिए बहुत अधिक उपयोग किया जाता है। एलडीयू कारक बनाने का प्रयास करें। या बेहतर, एल्गोरिदम के बारे में चिंता न करें, और इसके बजाय LAPACK का उपयोग करें। – Peter

+0

4 से कम चर के लिए, क्रैमर नियम नियम कोड इमो –

3

क्या आप एक ऐसे सॉफ्टवेयर पैकेज की तलाश में हैं जो काम करेगा या वास्तव में मैट्रिक्स संचालन कर रहा है और ऐसा और प्रत्येक चरण करता है?

पहला, मेरा एक सहकर्मी सिर्फ Ocaml GLPK का उपयोग करता था। यह GLPK के लिए सिर्फ एक रैपर है, लेकिन यह चीजों को सेट करने के कई चरणों को हटा देता है। ऐसा लगता है कि आपको सी में जीएलपीके के साथ रहना होगा, हालांकि। उत्तरार्द्ध के लिए, पुराने लेख को सहेजने के लिए स्वादिष्ट धन्यवाद, मैं थोड़ी देर पहले एलपी सीखने के लिए प्रयोग करता था, PDF। अगर आपको आगे बढ़ने में विशिष्ट सहायता की ज़रूरत है, तो हमें बताएं और मुझे यकीन है कि, मैं या कोई भी वापस भटक जाएगा और मदद करेगा, लेकिन, मुझे लगता है कि यह यहां से काफी आगे है। शुभ लाभ!

7

रैखिक समीकरणों की 3x3 प्रणाली के लिए मुझे लगता है कि यह आपके स्वयं के एल्गोरिदम को रोल करना ठीक होगा।

हालांकि, आपको सटीकता, शून्य या वास्तव में छोटी संख्याओं के विभाजन और असीमित कई समाधानों के बारे में चिंता करने की आवश्यकता हो सकती है। मेरा सुझाव मानक संख्यात्मक रैखिक बीजगणित पैकेज जैसे LAPACK के साथ जाना है।

1

व्यक्तिगत रूप से, मैं Numerical Recipes के एल्गोरिदम के आंशिक हूं। (मुझे सी ++ संस्करण का शौक है।)

यह पुस्तक आपको सिखाएगी कि एल्गोरिदम क्यों काम करते हैं, साथ ही आपको उन एल्गोरिदम के कुछ सुंदर-अच्छी तरह से डीबग किए गए कार्यान्वयन दिखाते हैं।

बेशक, आप अंधेरे से CLAPACK (मैंने इसे बड़ी सफलता के साथ उपयोग किया है) का उपयोग कर सकते हैं, लेकिन मैं पहले गॉसियन एलिमिनेशन एल्गोरिदम को हाथ से टाइप करूँगा ताकि कम से कम इस तरह के काम का एक बेहोश विचार हो इन एल्गोरिदम स्थिर बनाने में।

बाद में, यदि आप अधिक रोचक रैखिक बीजगणित कर रहे हैं, तो Octave के स्रोत कोड के चारों ओर देखकर बहुत सारे प्रश्नों का उत्तर दिया जाएगा।

3

Template Numerical Toolkit एनआईएसटी से ऐसा करने के लिए उपकरण हैं।

QR Decomposition का उपयोग करने के लिए अधिक विश्वसनीय तरीकों में से एक है।

यहां एक रैपर का एक उदाहरण दिया गया है ताकि मैं अपने कोड में "GetInverse (A, InvA)" को कॉल कर सकूं और यह इनवर्स को इनव में डाल देगा।

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

लाइब्रेरी में Array2D परिभाषित किया गया है।

+0

लिखने के लिए सबसे अच्छा है 'qr.solve (I)' में 'I' क्या है? – Ponkadoodle

2

आपके प्रश्न के शब्द से, ऐसा लगता है कि आपके पास अज्ञात से अधिक समीकरण हैं और आप विसंगतियों को कम करना चाहते हैं। यह आम तौर पर रैखिक प्रतिगमन के साथ किया जाता है, जो विसंगतियों के वर्गों के योग को कम करता है। डेटा के आकार के आधार पर, आप इसे स्प्रेडशीट में या सांख्यिकीय पैकेज में कर सकते हैं। आर एक उच्च गुणवत्ता वाला, मुफ्त पैकेज है जो कई अन्य चीजों के बीच रैखिक प्रतिगमन करता है। रैखिक प्रतिगमन (और बहुत सारे गॉचा) के लिए बहुत कुछ है, लेकिन सरल मामलों के लिए यह सरल है। यहां आपके डेटा का उपयोग करके एक आर उदाहरण है। ध्यान दें कि "टीएक्स" आपके मॉडल में अवरोध है।

> 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 
2

रन-टाइम दक्षता के मामले में, दूसरों की तुलना में बेहतर मैं जवाब दे दिया है आप हमेशा चर के रूप में समीकरण का एक ही नंबर होगा, तो मैं Cramer's rule तरह के रूप में इसे लागू करने में आसान है। बस एक मैट्रिक्स के निर्धारक की गणना करने के लिए एक फ़ंक्शन लिखें (या पहले से लिखे गए एक का उपयोग करें, मुझे यकीन है कि आप वहां से एक को ढूंढ सकते हैं), और दो मैट्रिस के निर्धारकों को विभाजित करें।

14

आप इसे एक प्रोग्राम के साथ ठीक उसी तरह हल कर सकते हैं जैसे आप इसे हाथ से हल करते हैं (गुणा और घटाव के साथ, फिर समीकरणों में परिणाम खिला रहे हैं)। यह सुंदर मानक माध्यमिक-विद्यालय स्तर के गणित है।

-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 = 0.0785250 
b = -0.1806125 
c = -41.6375875 

आप इन मूल्यों को ए, बी और सी में प्लग, तो आप वे सही कर रहे हैं मिल जाएगा।

चाल एक साधारण 4x3 मैट्रिक्स का उपयोग करना है जो 3x2 मैट्रिक्स में बदल जाती है, फिर एक 2x1 जो "ए = एन" है, एन वास्तविक संख्या है। एक बार आपके पास यह हो जाने के बाद, आप इसे अगले मैट्रिक्स में एक और मान प्राप्त करने के लिए फ़ीड करते हैं, फिर उन दो मानों को अगले मैट्रिक्स में तब तक फ़ीड करें जब तक आप सभी चर हल नहीं कर लेते।

बशर्ते आपके पास एन विशिष्ट समीकरण हों, तो आप हमेशा एन चर के लिए हल कर सकते हैं। मैं अलग कहते हैं कि इनमें नहीं हैं क्योंकि:

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

वे ही समीकरण दो से गुणा किया तो आप उन लोगों से कोई समाधान नहीं प्राप्त कर सकते हैं - सही लेकिन बेकार बयान के साथ दो तो घटाकर पत्ते आप से पहले गुणा :

0 = 0 + 0 

उदाहरण के द्वारा, यहाँ कुछ सी कोड है कि युगपत समीकरण है कि आप अपने प्रश्न में रखा रहे हैं बाहर काम करता है।सबसे पहले कुछ आवश्यक प्रकार, चर, एक समीकरण बाहर मुद्रण के लिए एक समर्थन समारोह, और 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; 

इसके बाद, दो अज्ञात के साथ दो समीकरणों को तीन अज्ञात के साथ तीन समीकरणों की कमी:

// 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 (""); 

इसके बाद, एक अज्ञात के साथ एक समीकरण को दो अज्ञात के साथ दो समीकरणों की कमी:

// 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 के साथ अज्ञात मान को आसानी से काम कर सकते हैं। फिर, एक बार जब आप उस मूल्य को समझ लेंगे, तो इसे दो अज्ञातों के साथ समीकरणों में से एक में बदलें और दूसरे मान को कार्य करें। फिर मूल समीकरण में से एक में उन दोनों (अब जाना जाता है) अज्ञात स्थानापन्न और अब आप सभी तीन अज्ञात के लिए मान हैं:

// 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) 
6

Microsoft Solver Foundation पर एक नज़र डालें।

इसके साथ

आप इस तरह कोड लिख सकते हैं:
=== सॉल्वर फाउंडेशन सेवा रिपोर्ट ===
Datetime: 04/20/2009 23

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

यहाँ उत्पादन होता है: 29:55
मॉडल का नाम: डिफ़ॉल्ट
क्षमताओं का अनुरोध: एल.पी.
का समाधान समय (मिलीसेकंड): 1027
कुल अवधि (मिलीसेकंड): 1414
सॉल्वर चयनित इष्टतम: Microsoft.SolverFoundation.Solvers.SimplexSolver
निर्देशों:
Microsoft.SolverFoundation.Services.Directive
एल्गोरिथ्म: प्राइमल
अंकगणित: हाइब्रिड
मूल्य निर्धारण (सटीक): चूक
1,363,210 समापन स्थिति का समाधान मूल्य निर्धारण (डबल): SteepestEdge
आधार: सुस्त
धुरी गणना: 3
=== समाधान विवरण ===
लक्ष्य:
+०१२३५०९५२९९३ निर्णय:
एक: .0785250000000004
ख: -0.180612500000001
c: -41,6375875

+0

तो, हम इससे किस संख्यात्मक स्थिरता गुणों की अपेक्षा कर सकते हैं? चूंकि यह खुला स्रोत नहीं है, इसलिए यह उचित परिश्रम के साथ आना चाहिए, एलएपीएक्स जैसे मुख्यधारा के पुस्तकालयों के खिलाफ बेंचमार्क, मालिकाना समाधान के साथ जाने के लिए कुछ महत्वपूर्ण लाभ होना चाहिए। –

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 
+0

तो क्या होगा यदि 'ए (एन, एन)' शून्य है? –