24小时热门版块排行榜    

查看: 3492  |  回复: 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的回帖

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的回帖
查看全部 18 个回答

东平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的回帖

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的回帖

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的回帖
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[博后之家] 售SCI一区文章,我:8O5.5.1.O5.4,科目全,可伽急 +3 TZqyNTBa07NP 2026-11-13 3/150 2026-11-15 10:19 by xSrQNSSjcgoX
[考研] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 SX5mJmPWDfOO 2026-11-13 3/150 2026-11-15 07:56 by xSrQNSSjcgoX
[论文投稿] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +4 SX5mJmPWDfOO 2026-11-13 4/200 2026-11-15 07:45 by xSrQNSSjcgoX
[考博] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +3 XJXMuRjD7cK7 2026-11-13 3/150 2026-11-15 07:23 by xSrQNSSjcgoX
[公派出国] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +3 uKOSCocCb6Fm 2026-11-13 3/150 2026-11-15 04:38 by xSrQNSSjcgoX
[考研] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +8 zSbvVX9JFl9J 2026-11-09 9/450 2026-11-13 11:11 by Equinoxhua
[硕博家园] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +5 C6CRn66RYYXS 2026-11-10 9/450 2026-11-13 11:02 by Equinoxhua
[基金申请] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +6 C6CRn66RYYXS 2026-11-10 7/350 2026-11-13 10:55 by Equinoxhua
[公派出国] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +4 cdqGD6ftHmbl 2026-11-11 5/250 2026-11-12 12:28 by ECLTXnJnRSXf
[公派出国] 售SCI一区T0P文章,我:8.O.55.1.O.54,科目齐全,可+急 +3 Z7qfCsdIiD1T 2026-11-10 7/350 2026-11-12 11:31 by HajghKOMynmz
[考博] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +3 Z7qfCsdIiD1T 2026-11-10 5/250 2026-11-12 11:25 by HajghKOMynmz
[考研] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +7 czAiolzR0IJC 2026-11-10 9/450 2026-11-12 11:00 by HajghKOMynmz
[教师之家] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +4 czAiolzR0IJC 2026-11-10 7/350 2026-11-12 11:00 by HajghKOMynmz
[考研] 售SCI文章,我:8O.5.5.1O.54,科目全,可十急 +7 0TEMgxdT1eIw 2026-11-09 9/450 2026-11-12 09:24 by HajghKOMynmz
[公派出国] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +4 Z7qfCsdIiD1T 2026-11-10 6/300 2026-11-12 07:09 by 9fOTqj074tIg
[公派出国] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +3 czAiolzR0IJC 2026-11-10 5/250 2026-11-12 06:36 by 9fOTqj074tIg
[教师之家] 售SCI一区T0P文章,我:8O.55.1.O.5.4,科目齐全,可+急 +5 FChKX8meu1ds 2026-11-09 10/500 2026-11-12 05:19 by 9fOTqj074tIg
[教师之家] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +4 FChKX8meu1ds 2026-11-09 7/350 2026-11-12 05:19 by 9fOTqj074tIg
[公派出国] 售SCI-T0P文章,我:8O.5.5.1.O.54,科目齐全,可+急 +5 zSbvVX9JFl9J 2026-11-09 8/400 2026-11-12 04:57 by 9fOTqj074tIg
[硕博家园] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +5 FChKX8meu1ds 2026-11-09 6/300 2026-11-11 17:03 by TgvkbZZl9tqB
信息提示
请填处理意见