征集求解线性方程组的数学软件

liangbch

安装了一个matlab 6.5,并尝试写了一个命令行脚本，现分享给大家。 用一个文本编辑器创建一个名为test1.m的文件，保存在C:\MATLAB6p5\work，文件内容如下：

1. A=[1,1,-5,-1,1;8,5,-1,0,4;3,-7,-5,-1,1;2,1,8,3,-1;0,6,2,7,-2];
2. b=[9;2;4;-8;15];
3. xx=A\b;
5. diary r.txt % 将所有输出保存至r.txt文件
6. diary on %打开diary开关
7. prot='A\b is ----------------';
9. disp('A=');
10. disp(A);
11. disp('b=');
12. disp(b);
13. disp(prot);
14. disp(xx);
15. diary off; %关闭diary开关
16. exit % 退出matlab

复制代码

在命令行运行matlab -r test1,则启动matlab,并求解方程组，保存结果到r.txt。

wayne

赞一个先， 不过，我得说我们一般都不是这样用MATLAB的~~ matlab -r test1最终还是会启动MATLAB图形界面的，到不如直接在的GUI界面里面的command windows里面输入三行代码来的爽快：

1. A=[1,1,-5,-1,1;8,5,-1,0,4;3,-7,-5,-1,1;2,1,8,3,-1;0,6,2,7,-2];
2. b=[9;2;4;-8;15];
3. xx=A\b;

复制代码

另外 MATLAB可以和c++进行混合编程的。。。

wayne

我的解法是消元法，将系数矩阵化为上三角矩阵。 liangbch 发表于 2010-3-1 16:31

一般的计算方法书籍都会讲到线性方程组求解的多种算法，消元法最基础。 我记得还有一个LU分解的方法

wayne

我来给一个gsl的版本：

1. #include <stdio.h>
2. #include <gsl/gsl_linalg.h>
3. int main ()
4. {
5. double a_data[] = {1,1,-5,-1,1,8,5,-1,0,4,3,-7,-5,-1,1,2,1,8,3,-1,0,6,2,7,-2};
6. double b_data[] = {9,2,4,-8,15};
8. gsl_matrix_view m = gsl_matrix_view_array (a_data, 5, 5);
9. gsl_vector_view b= gsl_vector_view_array (b_data, 5);
10. gsl_vector *x = gsl_vector_alloc (5);
11. int s;
12. gsl_permutation * p = gsl_permutation_alloc (5);
13. gsl_linalg_LU_decomp (&m.matrix, p, &s);
14. gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);
15. printf ("x = \n");
16. gsl_vector_fprintf (stdout, x, "%g");
18. gsl_permutation_free (p);
19. gsl_vector_free (x);
20. return 0;
21. }

复制代码

编译运行： H:\tkpng>gcc linear.c -Iinclude -Llib -lgsl -lgslcblas -o gsl H:\tkpng>gsl.exe x = 1.5024 1.0006 -2.30719 0.706587 -4.33234

liangbch

我成功地使用VC++6.0跑通了楼上的程序，下面列出我的方法 1)下载WinGsl-Lib-1.4.02.zip from http://www6.in.tum.de/~kiss/WinGsl.htm (2)解压到C:\，WinGsl中的库文件的扩展名是.lib,可用于Visual C++ 6.0 的编译器 (3)设置Visual C++ 6.0编译环境： 一、将 C:\WinGsl\Bin中的WinGsl.dll和WinGslD.dll复制到C:\Program Files\Microsoft Visual Studio\VC98\Bin下； 二、打开VC，tools－>options－>directories，将解压目录下的lib，gsl分别加入到库文件和头文件的搜索路径中， 即, 将C:\WINGSL加入到Include files 搜索路径，将C:\WINGSL\LIB加入到Library files. (4) 在新建的工程中将你所使用到的库文件加入到 Link->Object/library modules 中，我在debug模式下添加的是WinGslLib_md.lib，release 模式下添加的是WinGslLib_md.lib，不过好像用错了文件也可以工作。 参考自：DIY部落(http://www.diybl.com/course/3_pr ... 2008326/107387.html)

liangbch

还有一个问题，14#的数据类型是double型，不知能否使用GMP中的高精度数据类型。

liangbch

12# wayne 11楼的命令行用法 适用求解大量方程组，如果我们有10000个方程组，写一个批处理，采用11楼的的做法很快可以得到答案。如果用GUI则麻烦的多。

liangbch

　在我的应用中，有好多方程组是线性相关的。对于这样的方程组，在消元的过程中，某一个方程的所有系数会同时为0。这时，程序就可以确定方程组无解。 　我特意挑选了一个线性相关的方程组，使用 14# 的代码，用GLS求解，结果发现他依然给出一个解，当然这个解错误的。 　今天凌晨，检查了一下我的程序，在修复了一个bug后发现,整数溢出的概率非常低，再也没有看到整数溢出的报告。所以这个主题贴实际上对我的应用没有多大帮助。不过这个主题贴让我们学到很多，欢迎继续讨论，谢谢所有回复的人。

wayne

12# wayne 11楼的命令行用法 适用求解大量方程组，如果我们有10000个方程组，写一个批处理，采用11楼的的做法很快可以得到答案。如果用GUI则麻烦的多。 liangbch 发表于 2010-3-1 22:29

MATLAB是一个很强大的工具，即便是GUI也可以批处理的啊 我自己写了一个批处理的脚本，数据是随机产生的，有1000个方程组；你只需把代码保存到一个test.m文件里，在文件的焦点下按F5运行，或者在command windows里直接输入test也可以运行

1. %% 产生测试数据：1000个随机的线性方程组 Ax=b, A为5*5方阵,b为5*1列向量
2. for ii=1:1000
3. data(ii).A=randi([-2,8],5);
4. data(ii).b=randi([-2,8],5,1);
5. end
6. %% 批量解线性方程组,解统一存在结构表 data的x域里
7. for ii=1:1000
8. data(ii).x=data(ii).A\data(ii).b;
9. end
11. %% 把数据写入文件leqs.txt
12. file='leqs.txt'
13. for ii=1:10
14. dlmwrite(file, 'A=', 'delimiter', '','-append')
15. dlmwrite(file, data(ii).A, 'delimiter', '\t','-append')
16. dlmwrite(file, 'b=', 'delimiter', '','-append')
17. dlmwrite(file, data(ii).b,'delimiter', '\t', '-append')
18. dlmwrite(file, 'x=', 'delimiter', '','-append')
19. dlmwrite(file, data(ii).x, 'delimiter', '\t','-append')
20. dlmwrite(file, '================================','delimiter', '','-append')
21. end

复制代码

liangbch

19# wayne 你可能没有理解我的意思，我的问题是需要求解的方程组有好多个， 比如说有1000个，我已经用别的程序生成1000个m文件，如果你通过GUI的话，你需要 1.启动MATLAB 2.手动输入1000个命令来执行者1000个.m文件，如果我使用批处理方式来处理类似的问题，则可以创建一个名为 all.bat 批处理（内容如下），然后只需要在命令上输入all就可以了。 matlab -r 00001.m matlab -r 00002.m matlab -r 00003.m matlab -r 00004.m matlab -r 00005.m matlab -r 00006.m matlab -r 00006.m matlab -r 00007.m matlab -r 00008.m matlab -r 00009.m matlab -r 00010.m ...