Tôi cần lập trình giải quyết một hệ phương trình tuyến tính trong C, Mục tiêu C, hoặc (nếu cần) C++.Giải phương trình tuyến tính
Dưới đây là một ví dụ về các phương trình:
-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
Từ đó, tôi muốn để có được xấp xỉ tốt nhất cho a, b và tx.
Cramer's Rule và Gaussian Elimination hai thuật toán tốt, có mục đích chung (xem thêm Simultaneous Linear Equations). Nếu bạn đang tìm mã, hãy kiểm tra GiNaC, Maxima và SymbolicC++ (tùy thuộc vào yêu cầu cấp phép của bạn, tất nhiên).
EDIT: Tôi biết bạn đang làm việc ở đất C, nhưng tôi cũng phải đưa ra một từ tốt cho SymPy (một hệ thống đại số máy tính bằng Python). Bạn có thể học được rất nhiều từ các thuật toán của nó (nếu bạn có thể đọc một chút python). Ngoài ra, nó theo giấy phép BSD mới, trong khi hầu hết các gói toán học miễn phí là GPL.
thực sự, không phải quy tắc của cramer hay loại trừ gaussian đều rất tốt trong thế giới thực. không có thuộc tính số tốt, và không được sử dụng nhiều cho các ứng dụng nghiêm trọng. hãy thử hệ số LDU. hoặc tốt hơn, đừng lo lắng về thuật toán và sử dụng LAPACK thay thế. – Peter
cho các biến số nhỏ hơn 4, Quy tắc của Cramer là tốt nhất để viết mã giải mã imo –
Bạn đang tìm gói phần mềm sẽ thực hiện công việc hoặc thực sự thực hiện các thao tác ma trận và thực hiện các bước như vậy?
Người đầu tiên, đồng nghiệp của tôi vừa mới sử dụng Ocaml GLPK. Nó chỉ là một trình bao bọc cho GLPK, nhưng nó loại bỏ rất nhiều bước thiết lập mọi thứ. Có vẻ như bạn sẽ phải gắn bó với GLPK, trong C, mặc dù. Đối với thứ hai, nhờ ngon để tiết kiệm một bài báo cũ tôi đã từng học LP một thời gian trở lại, PDF. Nếu bạn cần trợ giúp cụ thể thiết lập thêm, hãy cho chúng tôi biết và tôi chắc chắn, tôi hoặc ai đó sẽ đi lang thang trở lại và giúp đỡ, nhưng, tôi nghĩ rằng nó khá thẳng về phía trước từ đây. Chúc may mắn!
Đối với hệ phương trình tuyến tính 3x3, tôi đoán bạn có thể triển khai các thuật toán của riêng mình.
Tuy nhiên, bạn có thể phải lo lắng về độ chính xác, chia cho số không hoặc số thực sự nhỏ và phải làm gì với vô số giải pháp. Đề nghị của tôi là đi với một gói đại số tuyến tính số chuẩn như LAPACK.
Cá nhân, tôi là một phần của thuật toán của Numerical Recipes. (Tôi thích phiên bản C++.)
Cuốn sách này sẽ cho bạn biết lý do tại sao các thuật toán hoạt động, đồng thời cho bạn thấy một số triển khai được gỡ lỗi khá tốt của các thuật toán đó. Tất nhiên, bạn chỉ có thể sử dụng một cách mù quáng CLAPACK (tôi đã sử dụng nó với thành công lớn), nhưng trước tiên tôi sẽ gõ một thuật toán loại bỏ Gaussian để ít nhất có một ý tưởng mờ nhạt về loại công việc đã biến mất làm cho các thuật toán này ổn định.
Sau đó, nếu bạn đang làm đại số tuyến tính thú vị hơn, nhìn xung quanh mã nguồn của Octave sẽ trả lời rất nhiều câu hỏi.
Template Numerical Toolkit từ NIST có các công cụ để thực hiện điều đó.
Một trong những cách đáng tin cậy hơn là sử dụng QR Decomposition.
Dưới đây là ví dụ về trình bao bọc để tôi có thể gọi "GetInverse (A, InvA)" trong mã của tôi và nó sẽ đưa ngược vào InvA.
void GetInverse(const Array2D<double>& A, Array2D<double>& invA)
{
QR<double> qr(A);
invA = qr.solve(I);
}
Mảng 2D được định nghĩa trong thư viện.
'I' trong' qr.solve (I) 'là gì? – Ponkadoodle
Từ từ ngữ của câu hỏi của bạn, có vẻ như bạn có nhiều phương trình hơn là ẩn số và bạn muốn giảm thiểu những mâu thuẫn. Điều này thường được thực hiện với hồi quy tuyến tính, giúp giảm thiểu tổng của các ô vuông không nhất quán. Tùy thuộc vào kích thước của dữ liệu, bạn có thể làm điều này trong một bảng tính hoặc trong một gói thống kê. R là một gói chất lượng cao, miễn phí có hồi quy tuyến tính, trong số rất nhiều thứ khác. Có rất nhiều để hồi quy tuyến tính (và rất nhiều của gotcha), nhưng vì nó đơn giản để làm cho các trường hợp đơn giản. Đây là một ví dụ R sử dụng dữ liệu của bạn. Lưu ý rằng "tx" là chặn cho mô hình của bạn.
> 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
Xét về hiệu quả thời gian chạy, những người khác đã trả lời tốt hơn so với I. Nếu bạn luôn luôn sẽ có cùng một số phương trình như biến, tôi thích Cramer's rule vì nó rất dễ dàng để thực hiện. Chỉ cần viết một hàm để tính toán yếu tố quyết định của một ma trận (hoặc sử dụng một hàm đã được viết, tôi chắc chắn bạn có thể tìm ra một ma trận), và chia các yếu tố quyết định của hai ma trận.
Bạn có thể giải quyết vấn đề này bằng chương trình giống hệt cách bạn giải quyết bằng tay (với phép nhân và phép trừ, sau đó cho kết quả đưa vào phương trình). Đây là toán học cấp trung học khá chuẩn.
-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)
Vì vậy, bạn kết thúc với:
a = 0.0785250
b = -0.1806125
c = -41.6375875
Nếu bạn cắm các giá trị trở lại vào A, B và C, bạn sẽ thấy họ là đúng.
Bí quyết là sử dụng ma trận 4x3 đơn giản làm giảm lần lượt thành ma trận 3x2, sau đó là 2x1 là "a = n", n là số thực. Một khi bạn có điều đó, bạn nạp nó vào ma trận tiếp theo để lấy một giá trị khác, sau đó hai giá trị đó vào ma trận tiếp theo cho đến khi bạn đã giải quyết tất cả các biến.
Với điều kiện bạn có N phương trình riêng biệt, bạn luôn có thể giải quyết các biến N. Tôi nói rõ ràng vì hai đây không phải là:
7a + 2b = 50
14a + 4b = 100
Họ là những phương trình cùng nhân với hai nên bạn không thể có được một giải pháp từ họ - nhân đầu tiên bởi hai sau đó trừ đi lá bạn với báo cáo kết quả đúng nhưng vô dụng :
0 = 0 + 0
Bằng một ví dụ, đây là một số mã C mà làm việc ra các phương trình đồng thời là bạn đang đặt trong câu hỏi của bạn.Đầu tiên một số loại cần thiết, các biến, một chức năng hỗ trợ cho in ra một phương trình, và khi bắt đầu 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;
Tiếp theo, việc giảm trong ba phương trình với ba ẩn số để hai phương trình với hai ẩn số:
// 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 ("");
Tiếp theo, việc giảm của hai phương trình với hai ẩn số với một phương trình với một không rõ:
// 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 ("");
Bây giờ chúng ta có một công thức của loại number1 = unknown * number2, chúng tôi chỉ có thể tìm ra giá trị không xác định với unknown <- number1/number2. Sau đó, khi bạn đã tìm ra giá trị đó, hãy thay thế nó thành một trong các phương trình với hai ẩn số và tính ra giá trị thứ hai. Sau đó, thay cả những ẩn số (nay là nổi tiếng) vào một trong các phương trình ban đầu và bây giờ bạn có giá trị cho tất cả ba ẩn số:
// 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;
}
Sản lượng mã mà phù hợp với tính toán trước đó trong câu trả lời này:
>: -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)
Hãy xem qua số Microsoft Solver Foundation.
Với nó, bạn có thể viết code như thế này:
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);
Đây là kết quả:
=== Solver Foundation Dịch vụ Báo cáo ===
Datetime: 2009/04/20 23: Tên 29:55
mẫu: Mặc định
Khả năng yêu cầu: LP
Giải quyết thời gian (ms): 1027
Tổng thời gian (ms): 1414
Giải quyết Tình trạng Hoàn thành: Optimal
Solver chọn: Microsoft.SolverFoundation.Solvers.SimplexSolver
Chỉ thị:
Microsoft.SolverFoundation.Services.Directive
Thuật toán: Primal
Arithmetic: Hybrid
giá (chính xác): Mặc định
giá (kép): SteepestEdge
Cơ sở: Slack
Pivot Count: 3
=== giải pháp chi tiết ===
Mục tiêu:
Quyết định:
a: 0,0785250000000004
b: -0,180612500000001
c: -41,6375875
Vậy, chúng ta có thể mong đợi những tính chất ổn định số nào từ điều này? Vì đây không phải là mã nguồn mở, nên nó phải đi kèm với sự tích cực, điểm chuẩn so với các thư viện chính như LAPACK, vv. Phải có một số lợi thế đáng kể để vượt qua việc phải đi với một giải pháp độc quyền. –
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
Vậy điều gì xảy ra nếu 'A (n, n)' là số không? –
Những người khác đã trả lời câu này, nhưng hãy kiểm tra các cuốn sách * Numerical Analysis: Toán học của khoa học máy tính * bởi Kincaid và Cheney. Cuốn sách chủ yếu là về việc giải các hệ phương trình khác nhau. – Matthew