| 查看: 2161 | 回复: 5 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
[求助]
病态矩阵求解
|
|||
|
病态线性方程组阵Ax=B，A为19*19的矩阵，条件数E11，行列式2.3E90。 （1）是不是这类病态矩阵直接采用A\B或者inv(A)*B结果偏差都比较大？试过用符号解再转换为数值和直接inv(A)*B，发现二者结果误差很小。对此问题比较疑惑。 （2）矩阵较大，不好贴出，A和B均在附件中，如果有其它高精度的求解方法（如tikhonov正则化方法，在网上下载了tikhonov源代码，采用后结果和A\B相比增大了1000倍，不知道是我自己写的调用代码有问题还是其它问题）求大神赐教。 [ Last edited by eudemonia333 on 2013-4-22 at 21:03 ] |
回复此楼» 本帖附件资源列表
-
欢迎监督和反馈：小木虫仅提供交流平台，不对该内容负责。
本内容由用户自主发布，如果其内容涉及到知识产权问题，其责任在于用户本人，如对版权有异议，请联系邮箱：xiaomuchong@tal.com - 附件 1 : A.xls
- 附件 2 : B.xls
2013-04-22 20:58:46, 16.5 K
2013-04-22 20:58:49, 13.5 K
» 收录本帖的淘帖专辑推荐
 精彩 |
» 猜你喜欢
070300化学279求调剂
已经有19人回复
-
材料考研调剂
已经有19人回复
-
272分材料子求调剂
已经有41人回复
-
求调剂，262机械专硕
已经有7人回复
-
266求调剂，一志愿哈工程电子信息，本科获多项国奖和省奖
已经有6人回复
-
0860004 求调剂 309分
已经有8人回复
-
270求调剂
已经有5人回复
-
347求调剂
已经有3人回复
-
280求调剂
已经有9人回复
-
药学305求调剂
已经有5人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 大家对病态矩阵是怎么理解的啊
已经有11人回复
- 跪求 在协方差矩阵非正定情况下，多元正态分布的概率分布函数求解方法！
已经有4人回复
- 求Matlab高手解决线性方程组的迭代求解问题
已经有17人回复
- 求大神帮忙解病态线性方程组
已经有5人回复
- 寻研究反问题的虫友？？？？
已经有22人回复
- 如何用matlab求解矩阵系数的二阶微分方程
已经有13人回复
- 向saladin983致敬，也请教问题
已经有6人回复
- 矩阵相乘，逆向求解，这种方法可行吗？
已经有8人回复
- inv函数
已经有3人回复
- 如何用c++来求解矩阵的相关问题。
已经有12人回复
- Catia机械有限元分析到最后出错提示
已经有7人回复
- 如何得到一个奇异矩阵的线性无关的行号
已经有3人回复
- 请教一下对某一个矩阵如何求解近似对称轴？
已经有3人回复
- 大型非对称稀疏矩阵求解
已经有9人回复
- 【求助】fortran求解矩阵
已经有7人回复
- 【求助】求解一道有关泡利矩阵的量子力学习题
已经有12人回复
- 【求助】Riccati传递矩阵法求解临界转速（有部分程序）
已经有5人回复
3楼2013-04-26 16:49:33
chunbendan
铜虫 (初入文坛)
- 应助: 2 (幼儿园)
- 金币: 237.9
- 红花: 1
- 帖子: 18
- 在线: 5.5小时
- 虫号: 2251057
- 注册: 2013-01-20
- 性别: GG
- 专业: 计算数学与科学工程计算
【答案】应助回帖
★
lovibond: 金币+1, 鼓励新虫应助 2013-04-27 12:12:00
lovibond: 金币+1, 鼓励新虫应助 2013-04-27 12:12:00
|
病态矩阵求解需要特别注意。matlab里直接用inv(A)*B肯定是不行的，误差超级大。用A\B不知道有没有内部进行处理，也不建议用。建议自己写算法，奇异值修正或者是迭代法都可以。我这里写两个简单的迭代法程序供楼主参考： Jacobi迭代的函数文件：XJGJacobi.m function [x1,n]=XJGJacobi(A,b,x0,eps) % [x1,n]=XJGJacobi(A,b,x0,eps)使用JACOBI迭代法进行迭代求解方程A*x=b % A，b 输入量，方程的矩阵和向量 % x0 输入量，迭代初始值 % eps 精度，比如0.0001 % x1 输出，结果 % n 输出，迭代次数 D=diag(diag(A)); B=D\(D-A); n=1; % 初始迭代 x1=B*x0+D\b; while norm(x1-x0,inf)>eps %终止条件 n=n+1; x0=x1; x1=B*x0+D\b; end GAUSS_SEIDIAL迭代的函数文件：XJGGauss.m function [x1,n]=XJGGauss(A,b,x0,eps) % [x1,n]=XJGGauss(A,b,x0,eps) % 本函数使用GAUSS-SEDIAL迭代法进行迭代求解线性方程A*x=b % A，b 方程中的矩阵和向量 % x0 初始迭代值 % eps 精度 % 2011年12月5日 18:47:56 LD=tril(A,0); U=-triu(A,1); B=LD\U; f=LD\b; x1=B*x0+f; n=1; %迭代次数 while norm(x1-x0,inf)>eps %终止条件 n=n+1; x0=x1; x1=B*x0+f; end return ％我随便写了A和B，你的矩阵太大了 %用于求解的脚本文件： A=[6,2,1,-2 2,5,0,-2 -2,0,8,5 1,3,2,7]; b=[4;7;-1;0]; x0=zeros(4,1); %初始值 eps=1e-5; %精度 % Jacobi迭代法求解并显示结果 [x1,n]=XJGJacobi(A,b,x0,eps); disp('Jacobi迭代的结果为：') disp(x1) disp('Jacobi迭代的次数为：') disp(n) % GAUSS-SEDIAL迭代法求解并显示结果 [x1,n]=XJGGauss(A,b,x0,eps); disp('Gauss-Sedial迭代的结果为：') disp(x1) disp('Guass-Sedial迭代次数为：') disp(n) |
2楼2013-04-26 09:27:39
chunbendan
铜虫 (初入文坛)
- 应助: 2 (幼儿园)
- 金币: 237.9
- 红花: 1
- 帖子: 18
- 在线: 5.5小时
- 虫号: 2251057
- 注册: 2013-01-20
- 性别: GG
- 专业: 计算数学与科学工程计算
4楼2013-04-26 23:24:09
chunbendan
铜虫 (初入文坛)
- 应助: 2 (幼儿园)
- 金币: 237.9
- 红花: 1
- 帖子: 18
- 在线: 5.5小时
- 虫号: 2251057
- 注册: 2013-01-20
- 性别: GG
- 专业: 计算数学与科学工程计算
【答案】应助回帖
|
上面两个程序我有地方偷懒了，直接用的\,不改了。我查了下，用A\B是可以的，其实这里的\运用的就是列主元消去法。 列主元素消去法是为控制舍入误差而提出来的一种算法,在Gauss消去法的消元过程中,若出现a=0,则消元无法进行,即使其不为0,但很小,把它作为除数,就会导致其他元素量级的巨大增长和舍入误差的扩散,最后使计算结果不可靠.而使用列主元素消去法计算,基本上能控制舍入误差的影响,并且选主元素比较方便。就是说列主元消去法能比较好的解决你问题里的误差。 我自己写了个列主元消去法的程序，算出来跟A\B的是一样的。误差在1e－10数量级，很精确了。 >> A\B ans = 1.0e+06 * -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 1.6746 1.5890 0.0015 >> A*ans-B ans = 1.0e-10 * 0.1792 0 0.0045 -0.0011 0.0057 0.0182 0.0001 0.0000 0.0364 0.1501 0 -0.0040 0.0169 0.0448 -0.0103 -0.0016 -0.0002 -0.0003 0.0002 function x=plufact(A,b) % x=plufact(A,b) 该函数用来进行线性方程组A*x=b的带主元的LU分解（高斯消去法） % Input A -N*N矩阵 % b -N*1向量 % output x -N*1向量，是方程组的解 % 肖金光编于2011年10月7日 09:31:12 [M,N]=size(A); if(M~=N) error('矩阵A的维数不匹配'); end A=cat(2,A,b); x=zeros(N,1); for k=1:N-1 % 找主元所在行 [max1,j] =max(abs(A(k:N,k))); if(max1==0) error('A是奇异矩阵，无唯一解'); end % 选定主元，交换行 M=A(k,  ;A(k, =A(j+k-1,;A(j+k-1, =M;clear M; % 进行消元 M=A(k+1:N,k)/A(k,k); A(k+1:N, =-M*A(k,+A(k+1:N,; end % U=A(:,1:N) x(N)=A(N,N+1)/A(N,N); for k=N-1:-1:1 x(k)=(A(k,end)-A(k,k+1:N)*x(k+1:N))/A(k,k); end 用这个算出来的 >> x=plufact(A,B) x = 1.0e+06 * -0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 1.6746 1.5890 0.0015 >> A*x-B ans = 1.0e-10 * 0.1064 0.0182 -0.0136 -0.0028 -0.0051 0 0.0000 0.0000 0.0728 0.0080 -0.0728 0.0182 0.0230 -0.0278 -0.0098 0.0016 -0.0001 -0.0002 -0.0009 结果是一样的 |
5楼2013-04-27 11:52:51
