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

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;
   
4. 
   
5. diary r.txt % 将所有输出保存至r.txt文件
   
6. diary on %打开diary开关
   
7. prot='A\b is ----------------';
   
8. 
   
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};
   
7.           
   
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");
    
17.      
    
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
   
10. 
    
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
...