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

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 仿真模拟 » MATLAB » 求Matlab高手解决线性方程组的迭代求解问题
181/2返回列表12下一页
查看: 2676  |  回复: 17
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

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 ]
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
自强不息，厚德载物
1楼 2012-11-24 16:07:22
已阅   回复此楼   关注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
赞一下(4人) 回复此楼
活着为了实现价值
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)

不过结果解还是一样的。
赞一下(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 就可以了
赞一下(1人) 回复此楼
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
赞一下(1人) 回复此楼
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 的主题更新
181/2返回列表12下一页
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭