注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

AlexYoung

做好大家都必须要做的事,发展自己感兴趣的事,结束圣人所不齿的事...

 
 
 

日志

 
 

高斯消元及主素消元法说明  

2010-10-26 16:57:52|  分类: 课程练习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |



1.上三角系数矩阵的线性方程组的解法

方式一:回带法的for循环表示

X(N)=b(N)/A(N,N);

for k=N-1:-1:1

    a=0;

    for i=k+1:N

        a=a+A(k,i)*X(i);

    end

    X(k)=(b(k)-a)/A(k,k);

end

方式二:回带法的矢量性表示

X(N)=b(N)/A(N,N);

for k=N-1:-1:1

    X(k)=(b(k)-sum(A(k,k+1:N)*X(k+1:N)))/A(k,k);

end

 

2.高斯消元模块:

方式一:

for k=1:N-1%消元次数

    D=A;%暂存上一步消元之后的第k列的结果,在下次消元过程中要用到。

    for j=k:N+1%为了处理第k行的系数

        A(k,j)=A(k,j)/D(k,k);

        for i=k+1:N%该循环是为了处理后续行的系数

            A(i,j)=A(i,j)-D(i,k)*A(k,j);

        end

    end

end

方式二:

for k=1:N-1%消元次数

    D=A(:,k);%暂存上一步消元之后的第k列的结果,在下次消元过程中要用到。

    for j=k:N+1%为了处理第k行的系数

        A(k,j)=A(k,j)/D(k);

        for i=k+1:N%该循环是为了处理后续行的系数

            A(i,j)=A(i,j)-D(i)*A(k,j);

        end

    end

end

方法三:

for k=1:N-1%消元次数

    for j=k:N+1%为了处理第k行的系数

        a=A(k,k);

        A(k,j)=A(k,j)/a;

        for i=k+1:N%该循环是为了处理后续行的系数

            D(i)=A(i,k);[A.Young1]  %暂存上一步消元之后的第k列的>=k行的结果,待用。

            A(i,j)=A(i,j)-D(i)*A(k,j);

        end

    end

end

方法四:

for k=1:N-1%消元次数

    D=A(:,k);%暂存上一步消元之后的第k列的结果,在下次消元过程中要用到。

    for j=k:N+1%为了处理第k行的系数

        A(k,j)=A(k,j)/D(k);

A(k+1:N,j)=A(k+1:N,j)-D(k+1:N)*A(k,j);

    end

end

 

注:模块中的D只是消元时用的暂存变量,为了尽量减少存储空间,以上三种方法依次进行了优化。方法一,直接将整个系数矩阵暂存在A中,即DN*N维的;方法二,只将第k列的数值暂存在D中,即DN*1维的;方法三,只将要用到的第k列的>=k行的值暂存在D中,减少了数据的复制次数,但是这里实现不了,因为经过j=k之后,A(i,k)=0了,所以应该在这个循环进行之前进行值的转存,若放在上一循环中则又会因为引用角标超出范围而引起语法错误。故总的来说最好还是用第二种消元方式,方法四只是在方法二的基础上将最里层的一个for循环改写为向量的形式,书写简便了,并没有实质性的算法变化。而且作为可读性来说方法二更好。

3.主素消元模块

for k=1:N-1

[y,p]=max(abs(B(k:N,k)));

C=B(k,:);

B(k,:)=B(p+k-1,:);

B(p+k-1,:)=C;%找主素,换位置

for j=k+1:N

        m=B(j,k)/ B(k,k);

        B(j,k:N+1)=B(j,k:N+1)-m*B(k,k:N+1);

    end

end

说明:此消元方式并没有像高斯消元那样要将主对角线消为1,所以是从左向右,再从上到下的顺序进行消元。


 [A.Young1]想法很好,但是这里是实现不了的。如果是将系数矩阵和值列向量分开的话应该是可以的。

  评论这张
 
阅读(303)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017