【数值分析matlab】

数值代数 顺序高斯消元法 matlab

**
基本思想： 通过初等变换将方程组转化为同解的三角形方程组（消元），再通过 回代求解该方程组的解。

因此，通过顺序消元法，利用matlab求方程组的解需要经过两个步骤：1.消元 ；2.回代

Stept1: 消元（自上而下）
把系数矩阵化为上三角矩阵要写几个循环？
首先是整个矩阵的迭代，每实现一个k，ak-1,k-1前的元素要为0（通俗来说就是第几次消元）；然后每行每列具体到每个元素aij 进行操作，需要三个循环。（其实，利用matlab矩阵优点，可以写成一个或者两个for循环）

刚开始 我只写了两个循环，发现只是其中的某行元素实现了消元，明明遍历到每个元素，整个消元过程并没有转起来，原因在哪里？缺少迭代过程，这和一个矩阵所有元素相加不一样的！

function X=gaosi(A,b)
%用高斯消元法求解方程组
%A表示方程组系数矩阵，b为右端列向量
%x为解向
%% 例子
A=[1 1 1;-1 3 1;2 -6 1];  b=[6;4;-5];
%% 消元
n=length(A);a=[A,b];
for k=2:n
   for i=k:n
   a(i,k-1:n+1)=a(i,k-1:n+1)-a(k-1,k-1:n+1).*(a(i,k-1)/a(k-1,k-1)) ; % x=a;
   end
end

Stept2: 回代（自下而上）
这个就简单多了，先求解出来xn，然后倒序求解xj。只需要将求解过程分解，求和与除。

%% 回代
X=zeros(n,1); %写一个矩阵以便放入解向量；
X(n,1)=a(n,n+1)/a(n,n); %h放在这里的话会出错，要放在第一个循环后面，不然变成所有h的和
%此下求解分为两部分，需要写入一个求和
for i=(n-1):-1:1
    h=0;
    for j=i+1:n
      h=a(i,j)*X(j,1)+h;
      h;
    end
   X(i,1)=(a(i,n+1)-h)/a(i,i);
   X
end

好了，以上输入方程组就可以求解。但是要求ak,k不能等于0，而且还可能出现“大数吃小数”现象，要想避免这种现象，那就用到列主元高斯消元法。

附上，只用消元是只用一个循环的matlab代码（不得不感叹，通过矩阵写真方便）（此为书中代码）

function x=nagauss(a,b,flag)
%H1行 用途顺序高斯消去法解线性方程组ax=b;
%a:系数矩阵；b:右端列向量
%flag:若为0，则显示中间过程，否则不显示，默认值为0
%0：解向量
if nargin<0
    flag=0;
end
n=length(b);a=[a,b]
%消元
for k=1:(n-1)
    a((k+1):n,(k+1):(n+1))= a((k+1):n,(k+1):(n+1))-a((k+1):n,k)/a(k,k)*a(k,(k+1):(n+1));
    a(k+1:n,k)=zeros(n-k,1);
    if flag==0,a,end
end
%回代
x=zeros(n,1);
x(n)=a(n,n+1)/a(n,n);
for k=n-1:-1:1
    x(k)=(a(k,n+1)-a(k,(k+1):n)*x((k+1:n)))/a(k,k);
end

参考书：《数值分析及其MATLAB实验》姜健飞