应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 专业学科区 » 数学 » 计算数学 » 病态矩阵求解
61/1返回列表
查看: 2145  |  回复: 5
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

eudemonia333

新虫 (初入文坛)

[求助] 病态矩阵求解

病态线性方程组阵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
    • 2013-04-22 20:58:46, 16.5 K
  • 附件 2 : B.xls
    • 2013-04-22 20:58:49, 13.5 K

» 收录本帖的淘帖专辑推荐

精彩

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼 2013-04-22 20:59:04
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

chunbendan

铜虫 (初入文坛)

【答案】应助回帖

★
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
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

eudemonia333

新虫 (初入文坛)
引用回帖:
2楼: Originally posted by chunbendan at 2013-04-26 09:27:39
病态矩阵求解需要特别注意。matlab里直接用inv(A)*B肯定是不行的，误差超级大。用A\B不知道有没有内部进行处理，也不建议用。建议自己写算法，奇异值修正或者是迭代法都可以。我这里写两个简单的迭代法程序供楼主参 ...

谢谢！
上面两个程序不太适应我的矩阵，计算完出现了NaN，Warning: Matrix is singular to working precision. 直接采用inv没有上面两个错误。
赞一下 回复此楼
3楼2013-04-26 16:49:33
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

chunbendan

铜虫 (初入文坛)

【答案】应助回帖

引用回帖:
3楼: Originally posted by eudemonia333 at 2013-04-26 16:49:33
谢谢！
上面两个程序不太适应我的矩阵，计算完出现了NaN，Warning: Matrix is singular to working precision. 直接采用inv没有上面两个错误。...

我没用你的数据测试，应该修改一下能用，或者用列主元消去法也可以。我明天帮你弄一下
用inv处理奇异矩阵不行的

[ 发自手机版 http://muchong.com/3g ]
赞一下 回复此楼
4楼2013-04-26 23:24:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

chunbendan

铜虫 (初入文坛)

【答案】应助回帖

上面两个程序我有地方偷懒了，直接用的\,不改了。我查了下，用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
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

chunbendan

铜虫 (初入文坛)

【答案】应助回帖

寒，怎么变成表情了哦。 就是一个：加一个）
赞一下 回复此楼
6楼2013-04-27 11:55:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 eudemonia333 的主题更新
61/1返回列表
普通表情 龙 兔 虎 猫
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 303求调剂 +5 睿08 2026-03-17 7/350 2026-03-21 03:11 by JourneyLucky
[考研] 初始318分求调剂（有工作经验） +3 1911236844 2026-03-17 3/150 2026-03-21 02:33 by JourneyLucky
[考研] 279分求调剂 一志愿211 +11 chaojifeixia 2026-03-19 12/600 2026-03-21 01:49 by 星空星月
[考研] 华东师范大学-071000生物学-293分-求调剂 +3 研究生何瑶明 2026-03-18 3/150 2026-03-21 01:30 by JourneyLucky
[考研] 一志愿华南师大 070300（化学）304分求调剂 +3 0703武芊慧雪304 2026-03-18 3/150 2026-03-21 00:48 by JourneyLucky
[考研] A区线材料学调剂 +5 周周无极 2026-03-20 5/250 2026-03-20 21:33 by laoshidan
[考研] 求调剂 +3 @taotao 2026-03-20 3/150 2026-03-20 19:35 by JourneyLucky
[考研] 086500 325 求调剂 +3 领带小熊 2026-03-19 3/150 2026-03-20 18:38 by 尽舜尧1
[考研] 081700化工学硕调剂 +3 【1】 2026-03-16 3/150 2026-03-19 23:40 by edmund7
[考研] 材料，纺织，生物（0856、0710），化学招生啦 +3 Eember. 2026-03-17 9/450 2026-03-18 10:28 by Eember.
[考研] 0703化学336分求调剂 +6 zbzihdhd 2026-03-15 7/350 2026-03-18 09:53 by zhukairuo
[考研] 334求调剂 +3 志存高远意在机� 2026-03-16 3/150 2026-03-18 08:34 by lm4875102
[考博] 26博士申请 +3 1042136743 2026-03-17 3/150 2026-03-17 23:30 by 轻松不少随
[考研] 326求调剂 +5 上岸的小葡 2026-03-15 6/300 2026-03-17 17:26 by ruiyingmiao
[考研] 一志愿苏州大学材料工程（085601）专硕有科研经历三项国奖两个实用型专利一项省级立项 +6 大火山小火山 2026-03-16 8/400 2026-03-17 15:05 by 无懈可击111
[考研] 材料工程专硕274一志愿211求调剂 +6 薛云鹏 2026-03-15 6/300 2026-03-17 11:05 by 学员h26Tkc
[考研] [导师推荐]西南科技大学国防/材料导师推荐 +3 尖角小荷 2026-03-16 6/300 2026-03-16 23:21 by 尖角小荷
[考研] 0854控制工程 359求调剂 可跨专业 +3 626776879 2026-03-14 9/450 2026-03-16 17:42 by 626776879
[考研] 304求调剂 +4 ahbd 2026-03-14 4/200 2026-03-16 16:48 by 我的船我的海
[考研] 085601材料工程315分求调剂 +3 yang_0104 2026-03-15 3/150 2026-03-15 10:58 by peike
信息提示
关闭
请填处理意见
关闭