24小时热门版块排行榜    

查看: 2823  |  回复: 17

wangjunyahit

银虫 (初入文坛)

[求助] 求Matlab高手解决线性方程组的迭代求解问题

遇到一个线性方程组Ax=b的求解问题,使用Matlab直接求解但不确定结果是否正确,
想使用迭代解法求解验证结果,但苦于找不到合适的迭代解法和程序,希望各位大神
给予帮助,感激涕零!!!
方程组是这样的:
A=[-1.2e-21,1.4e6,1.4e6,-1.4e6,1.3e6;
    8.9e-22,-1.4e6,4.9e3,3.8e3,0;
    1.8e-22,1.4e3,-1.4e6,4.3e3,0;
    4.4e-23,3.2e2,1.2e3,-5.1e6,4.3e3;
    5.3e-23,0,0,1.3e3,-2.3e6]
b=[0;-1.8e-4;-1.3e-2;-1.5e-2;-4.0e-3]
直接求解如下:
x=A\b;
x=[  1.029393177754837e+20
     6.565808888911404e-08
     2.259822851544940e-08
     3.842185487752932e-09
     4.113382036189007e-09]

rcond(A)=1.694258513618419e-29;矩阵A可逆的条件数估值接近0,矩阵A条件不好,得到的解有误差
希望各位大神
给予帮助,感激涕零!!!

~EZ9ND[(A@F3IUUY9EJVQ4C.jpg

[ Last edited by wangjunyahit on 2012-11-24 at 16:10 ]
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

自强不息,厚德载物
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

fabaoyi

新虫 (初入文坛)

【答案】应助回帖

★ ★
感谢参与,应助指数 +1
xiegangmai: 金币+2, 谢谢参与 2012-11-25 19:58:30
你学过数值分析没有,这个可以用高斯迭代法和雅可比迭代法求解,不过虽然用迭代法,跑程序时还是有警告,可能求解还是不是非常可靠,我求解了一个,下面我把高斯迭代法那个我求解用的程序和运行结果贴在下面:
function gsdd(A,b,x0,fs,wucha,ddcs)
D=diag(diag(A));
U=-triu(A,1);
L=-tril(A,-1);
ljzx=inv(A)*b';
if det(D)==0
    disp('对角矩阵D奇异,迭代不能实现');
    return
end
B=(D-L)\U;f=(D-L)\b';
for k=1:ddcs
    x1=B*x0'+f;
    epsilon=norm(x1-x0',fs);
    wcx=norm(x1-x0',fs)/(norm(x1,fs)+eps);
    x0=x1';
    if (epsilon         break
    end
end
if (epsilon>wucha)||(wcx>wucha)||k==ddcs
   disp('   超过给定最大迭代次数,可能不收敛 ')
    return
end
disp('迭代次数:');
k
disp('迭代结果:');
x=x0'
disp('逆矩阵求解结果:');
x=ljzx
end
>> A=[-1.2e-21,1.4e6,1.4e6,-1.4e6,1.3e6;
     8.9e-22,-1.4e6,4.9e3,3.8e3,0;
     1.8e-22,1.4e3,-1.4e6,4.3e3,0;
     4.4e-23,3.2e2,1.2e3,-5.1e6,4.3e3;
     5.3e-23,0,0,1.3e3,-2.3e6];
>> b=[0;-1.8e-4;-1.3e-2;-1.5e-2;-4.0e-3]';
>> x0=[1 1 1 1 1];
>> gsdd(A,b,x0,1,10^(-15),1000)
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.
RCOND = 1.694259e-029.
> In gsdd at 5
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.
RCOND = 2.352941e-028.
> In gsdd at 10
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.
RCOND = 2.352941e-028.
> In gsdd at 10
迭代次数:

k =

   527

迭代结果:

x =

    1.029393177754839e+020
    6.565808888911417e-008
    2.259822851544943e-008
    3.842185487752933e-009
    4.113382036189012e-009

逆矩阵求解结果:

x =

    1.029393177754837e+020
    6.565808888911401e-008
    2.259822851544940e-008
    3.842185487752931e-009
    4.113382036189006e-009
活着为了实现价值
3楼2012-11-25 01:28:13
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

东平Modelon

金虫 (小有名气)

【答案】应助回帖

★ ★ ★ ★ ★
感谢参与,应助指数 +1
xiegangmai: 金币+2, 谢谢应助 2012-11-24 20:54:11
wangjunyahit: 金币+3, ★★★很有帮助 2012-11-25 11:12:49
矩阵的rcond显然可以改好啊。例如将第一列乘以4e27。

>> Anew = A; Anew(:,1) = A(:,1) *4e27;rcond(Anew)
    0.0103
>>x_1= Anew\b; x_1_real = x_1; x_1_real(1) = 4e27 * x_1(1)

不过结果解还是一样的。
2楼2012-11-24 20:35:20
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

change0618

铁杆木虫 (著名写手)

方丈大师

【答案】应助回帖


xiegangmai: 金币+1, 谢谢参与! 2013-04-27 23:43:15
对于刚性的矩阵,采用x=pinv(A)*b 就可以了
15楼2013-04-27 15:35:47
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

change0618

铁杆木虫 (著名写手)

方丈大师

【答案】应助回帖

★ ★ ★ ★
xiegangmai: 金币+2, 谢谢参与! 2013-04-27 23:43:35
wangjunyahit: 金币+2, ★★★很有帮助 2013-04-30 16:19:30
CODE:
>> A=[-1.070906744447462E-017   244564571.598802        303866255.019473      303856255.019473        303756255.019473  ;  
    9.735515836938967E-017   -437040200.438244        546168396.621813        415770404.263350      0.000000000000000E+000 ;
    1.947103215071310E-017   157778358.747585     -986775896.684545        473345946.492110       0.000000000000000E+000;
    4.867758037678274E-018   34697270.0918581        136741245.043259       -1333413850.81819        473345946.492110   ;   
     5.841309406796348E-018  0.000000000000000E+000  0.000000000000000E+000   136741245.043259      -778102201.511583 ];
b=[ 104348331231.704       -2072871.66710357      -154296552.906917       -178403143.022063       -48501206.8009638 ]';
>> x=pinv(A)*b

x =

   1.0e+02 *

  -0.000000000000000
   1.682414884915961
   0.668383858892869
   0.388902999746503
   0.336174005405975

>> vpa(x,30)

ans =

-0.0000000000000000000000285841068286876334276912171331
                         168.241488491596101084724068642
                         66.8383858892869398005132097751
                         38.8902999746503041933465283364
                         33.6174005405974725135820335709

17楼2013-04-27 15:48:18
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通回帖

wangjunyahit

银虫 (初入文坛)

引用回帖:
2楼: Originally posted by 东平Modelon at 2012-11-24 20:35:20
矩阵的rcond显然可以改好啊。例如将第一列乘以4e27。

>> Anew = A; Anew(:,1) = A(:,1) *4e27;rcond(Anew)
    0.0103
>>x_1= Anew\b; x_1_real = x_1; x_1_real(1) = 4e27 * x_1(1)

不过结果解

对的,我也这样试过,谢谢了
自强不息,厚德载物
4楼2012-11-25 09:28:56
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wangjunyahit

银虫 (初入文坛)

引用回帖:
2楼: Originally posted by 东平Modelon at 2012-11-24 20:35:20
矩阵的rcond显然可以改好啊。例如将第一列乘以4e27。

>> Anew = A; Anew(:,1) = A(:,1) *4e27;rcond(Anew)
    0.0103
>>x_1= Anew\b; x_1_real = x_1; x_1_real(1) = 4e27 * x_1(1)

不过结果解

我想问下将rcond改好之后是不是意味着方程的刚性减弱,不再是病态矩阵,这样的直接解法可以认为是正确的?
自强不息,厚德载物
5楼2012-11-25 09:33:37
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wangjunyahit

银虫 (初入文坛)

引用回帖:
3楼: Originally posted by fabaoyi at 2012-11-25 01:28:13
你学过数值分析没有,这个可以用高斯迭代法和雅可比迭代法求解,不过虽然用迭代法,跑程序时还是有警告,可能求解还是不是非常可靠,我求解了一个,下面我把高斯迭代法那个我求解用的程序和运行结果贴在下面:
function g ...

谢谢您的解答,非常感激!
我之前也用过雅克比迭代和高斯赛德尔迭代进行过求解另外一个类似的方程,但是没有收敛,不知道是不是初值选择的问题,麻烦您帮忙看下!谢谢!方程如下:
A=[-1.070906744447462E-017   244564571.598802        303866255.019473      303856255.019473        303756255.019473  ;  
   9.735515836938967E-017   -437040200.438244        546168396.621813        415770404.263350      0.000000000000000E+000 ;
   1.947103215071310E-017   157778358.747585     -986775896.684545        473345946.492110       0.000000000000000E+000;
   4.867758037678274E-018   34697270.0918581        136741245.043259       -1333413850.81819        473345946.492110   ;   
    5.841309406796348E-018  0.000000000000000E+000  0.000000000000000E+000   136741245.043259      -778102201.511583 ];
b=[ 104348331231.704       -2072871.66710357      -154296552.906917       -178403143.022063       -48501206.8009638 ]';
自强不息,厚德载物
6楼2012-11-25 11:12:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

fabaoyi

新虫 (初入文坛)

【答案】应助回帖

★ ★
xiegangmai: 金币+2, 帖中代码部分可使用“插入代码”功能 2012-11-25 19:58:47
引用回帖:
6楼: Originally posted by wangjunyahit at 2012-11-25 11:12:12
谢谢您的解答,非常感激!
我之前也用过雅克比迭代和高斯赛德尔迭代进行过求解另外一个类似的方程,但是没有收敛,不知道是不是初值选择的问题,麻烦您帮忙看下!谢谢!方程如下:
A=;
b=';...

这个用迭代法的话验证的话,没有哪一种迭代法是通用的,验证的话,你就用列主元消去法吧,一定程度上能控制舍入误差.下面且用列主元消去法求解,过程贴下面:
function gauss(A,b,delta)
%解线性方程组
B=[A b'];
[B,n]=disha(B);
X=A\b';
for k=1:n
    if abs(B(k,k))         disp('注意:detA=0!');
        return
    end
end
x=(B(:,1:n))\B(:,n+1);
disp('方程有有唯一解其解为:');
disp('i               x(i)');
for i=1:n
    disp(sprintf('%d          %.17e',i,x(i)));
end
    function [A hs]=disha(A)
        [hs,ls]=size(A);
        for i=1:hs-1
            s=A(:,i);
            [m,row]=max(abs(s(i:hs)));
            if A(i,i)~=m
                s1=A(i,;
                A(i,=A(row+i-1,;
                A(row+i-1,=s1;
            end
            for j=i:hs-1
                %   if A(j+1,i)==0
                %     continue
                %  end
                % r=A(i,i)/A(j+1,i);
                %A(j+1,=A(j+1,*r-A(i,;
                r=A(j+1,i)/A(i,i);
                A(j+1,=A(i,*r-A(j+1,;
            end
        end
    end
disp('逆矩阵求解结果: ');
for k=1:1:5
disp(sprintf('%d          %.17e',k,X(k)'));
end
end
>> gauss(A,b,10^(-18))
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 4.382649e-026.
> In gauss at 5
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 6.691293e-026.
> In gauss at 12
方程有有唯一解其解为:
i               x(i)
1          8.90763948846464080e+026
2          3.22058008969693620e+002
3          8.07073018970861680e+001
4          2.39314314743592430e+001
5          1.09550426143824530e+001
逆矩阵求解结果:
1          8.90763948846463800e+026
2          3.22058008969693560e+002
3          8.07073018970861540e+001
4          2.39314314743592430e+001
5          1.09550426143824530e+001
>>
活着为了实现价值
7楼2012-11-25 13:04:02
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

wangjunyahit

银虫 (初入文坛)

引用回帖:
7楼: Originally posted by fabaoyi at 2012-11-25 13:04:02
这个用迭代法的话验证的话,没有哪一种迭代法是通用的,验证的话,你就用列主元消去法吧,一定程度上能控制舍入误差.下面且用列主元消去法求解,过程贴下面:
function gauss(A,b,delta)
%解线性方程组
B=;
=disha(B ...

看不懂啊,能不能给下面的语句加点注释
[B,n]=disha(B);                %??
x=(B(:,1:n))\B(:,n+1)       %??
s1=A(i,;                         %??
                A(i,=A(row+i-1,;
                A(row+i-1,=s1;
我的matlab水平一般,多谢了,呵呵
自强不息,厚德载物
8楼2012-11-25 17:04:02
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

东平Modelon

金虫 (小有名气)


xiegangmai: 金币+1 2012-11-25 19:59:01
引用回帖:
5楼: Originally posted by wangjunyahit at 2012-11-25 02:33:37
我想问下将rcond改好之后是不是意味着方程的刚性减弱,不再是病态矩阵,这样的直接解法可以认为是正确的?...

显然如此。rcond (matlab帮助建议用cond(a))就是确定是否病态啊。
9楼2012-11-25 17:38:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

fabaoyi

新虫 (初入文坛)

【答案】应助回帖


xiegangmai: 金币+1, 谢谢参与 2012-11-26 19:08:26
引用回帖:
8楼: Originally posted by wangjunyahit at 2012-11-25 17:04:02
看不懂啊,能不能给下面的语句加点注释
=disha(B);                %??
x=(B(:,1:n))\B(:,n+1)       %??
s1=A(i,;                         %??
                A(i,=A(row+i-1,;
                A(row+i ...

[m,row]=max(abs(s(i:hs)));%找最大列主元,m为最大列主元,row为所在行
            if A(i,i)~=m%如果对角元不等于最大列主元
                s1=A(i,;%就将最大列主元与对角线上的元交换
                A(i,=A(row+i-1,;
                A(row+i-1,=s1;
            end
活着为了实现价值
10楼2012-11-25 21:55:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 wangjunyahit 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 279分求调剂 一志愿211 +9 chaojifeixia 2026-03-19 10/500 2026-03-20 12:29 by lature00
[考研] 招收调剂硕士 +4 lidianxing 2026-03-19 12/600 2026-03-20 12:25 by lidianxing
[考博] 招收博士1-2人 +3 QGZDSYS 2026-03-18 3/150 2026-03-20 11:58 by 呱呱呱呱叫
[考研] 085700资源与环境308求调剂 +10 墨墨漠 2026-03-18 11/550 2026-03-20 11:57 by 学员8dgXkO
[考研] 317求调剂 +4 申子申申 2026-03-19 8/400 2026-03-20 11:20 by 申子申申
[考研] 一志愿吉林大学材料学硕321求调剂 +11 Ymlll 2026-03-18 14/700 2026-03-20 10:46 by EBSD
[考研] 0703化学调剂 +10 妮妮ninicgb 2026-03-15 14/700 2026-03-19 22:59 by 学员8dgXkO
[考研] 一志愿北京化工大学0703化学318分,有科研经历,求调剂 +3 一瓶苯甲酸 2026-03-14 3/150 2026-03-19 15:17 by 尽舜尧1
[考研] 材料与化工求调剂 +7 为学666 2026-03-16 7/350 2026-03-19 14:48 by 尽舜尧1
[考研] 0703化学调剂 +5 pupcoco 2026-03-17 8/400 2026-03-19 13:58 by houyaoxu
[考研] 328求调剂,英语六级551,有科研经历 +4 生物工程调剂 2026-03-16 12/600 2026-03-19 11:10 by 生物工程调剂
[考研] 085601专硕,总分342求调剂,地区不限 +5 share_joy 2026-03-16 5/250 2026-03-18 14:48 by haxia
[考研] 297求调剂 +8 戏精丹丹丹 2026-03-17 8/400 2026-03-18 14:30 by laoshidan
[考博] 26申博 +4 八6八68 2026-03-16 4/200 2026-03-17 13:00 by 轻松不少随
[考研] 一志愿,福州大学材料专硕339分求调剂 +3 木子momo青争 2026-03-15 3/150 2026-03-17 07:52 by laoshidan
[考研] 304求调剂 +4 ahbd 2026-03-14 4/200 2026-03-16 16:48 by 我的船我的海
[考研] 321求调剂 +5 大米饭! 2026-03-15 5/250 2026-03-16 16:33 by houyaoxu
[考研] 277材料科学与工程080500求调剂 +3 自由煎饼果子 2026-03-16 3/150 2026-03-16 14:10 by 运气yunqi
[考研] 本科南京大学一志愿川大药学327 +3 麦田耕者 2026-03-14 3/150 2026-03-14 20:04 by 外星文明
[考研] 304求调剂 +7 7712b 2026-03-13 7/350 2026-03-13 21:42 by peike
信息提示
请填处理意见