1，插值算法（3种）：
（1）最邻近插值（近邻取样法）：
最邻近插值的的思想很简单，就是把这个非整数坐标作一个四舍五入，取最近的整数点坐标处的点的颜色。可见，最邻近插值简单且直观，速度也最快，但得到的图像质量不高。
最邻近插值法的MATLAB源代码为：
A = imread（'F:lena.jpg'）；%读取图像信息 imshow(A)；%显示原图 title（'原图128*128'）；
Row = size(A,1)；Col = size(A,2)；%图像行数和列数 nn=8;%放大倍数
m = round(nn*Row)；%求出变换后的坐标的最大值 n = round(nn*Col)；
B = zeros(m,n,3)；%定义变换后的图像
for i = 1 : m
for j = 1 : n
x = round(i/nn)；y = round(j/nn)；%最小临近法对图像进行插值
if x==0 x = 1;end
if y==0 y = 1;end
if x>Row x = Row;end
if y>Col y = Col;end B（i,j,:）= A（x,y,:）；
end end
B = uint8(B)；%将矩阵转换成8位无符号整数 figure;imshow(B)；
title（'最邻近插值法放大8倍1024*1024'）；
运行程序后，原图如图1所示：
图1
用最邻近插值法放大4倍后的图如图2所示：
图2
（2）双线性内插值法：
在双线性内插值法中，对于一个目的像素，设置坐标通过反向变换得到的浮点坐标为(i+u,j+v)，其中i、j均为非负整数，u、v为[0,1)区间的浮点数，则这个像素得值 f(i+u,j+v)可由原图像中坐标为(i,j)、(i+1,j)、(i,j+1)、(i+1,j+1)所对应的周围四个像素的值决定，即：
f(i+u,j+v)=(1-u)(1-v)f(i,j)+(1-u)vf(i,j+1)+ u(1-v)f(i+1,j)+ uvf(i+1,j+1)其中f(i,j)表示源图像(i,j)处的的像素值，以此类推。
这就是双线性内插值法。双线性内插值法计算量大，但缩放后图像质量高，不会出现像素值不连续的的情况。由于双线性插值具有低通滤波器的性质，使高频分量受损，所以可能会使图像轮廓在一定程度上变得模糊。
在MATLAB中，可用其自带的函数imresize()来实现双线性内插值算法。
双线性内插值算法的MATLAB源代码为：
A=imread（'F:lena.jpg'）；imshow(A)；
title（'原图128*128'）；
C=imresize（A,8,'bilinear'）；%双线性插值 figure;imshow(C)；
title（'双线性内插值法放大8倍1024*1024'）；
程序运行后，原图如图3所示：
图3
双线性内插值法放大8倍后的图如图4所示：
图4
（3）双三次插值法：
双三次插值法能够在很大程度上克服以上两种算法的不足，该算法计算精度高，但计算量大，它考虑一个浮点坐标(i+u,j+v)周围的16个邻点。
目的像素值f（i+u,j+v)可由如下插值公式得到：f(i+u,j+v)= [A] * [B] * [C] 其中[A]=[ S(u + 1)S(u + 0)S(u2)]；[C]=[ S(v + 1)S(v + 0)S(v2)]；而[B]是周围16个邻点组成的4*4的矩阵；S(x)是对 Sin(x*π）/x 的逼近。
在MATLAB中，可用其自带的函数imresize()来实现双三次插值算法。MATLAB源代码为：
A=imread（'F:lena.jpg'）；%读取原图像
D=imresize（A,8,'bicubic'）；%双三次插值放大8倍 figure;
T
imshow(D)；title（'三次卷积法放大8倍1024*1024'）；
MATLAB自带双三次插值法运行结果如图5所示：
图5
也可以自己编写双三次插值算法MATLAB代码如下：
clc,clear;
ff=imread（'F:lena.jpg'）；%读取图像到ff
k=8;%设置放大倍数 [m,n,color]=size(ff)；
f=zeros(m,n)；%将彩色图像ff转换为黑白图像f for i=1:m
for j=1:n
f(i,j)=ff(i,j)；
end end
a=f（1,:）；c=f（m,:）；%将待插值图像矩阵前后各扩展两行两列，共扩展四行四列 b=[f(1,1)，f(1,1)，f（:，1）'，f(m,1)，f(m,1)]；d=[f(1,n)，f(1,n)，f（:，n）'，f(m,n)，f(m,n)]；a1=[a;a;f;c;c]；a1'；
b1=[b;b;a1'；d;d]；f=b1'；f1=double(f)；
for i=1:k*m %利用双三次插值公式对新图象所有像素赋值 u=rem(i,k)/k;i1=floor(i/k)+2;A=[sw(1+u)sw(u)sw(1-u)sw(2-u)]；
for j=1:k*n
v=rem(j,k)/k;j1=floor(j/k)+2;C=[sw(1+v)；sw(v)；sw(1-v)；sw(2-v)]；
B=[f1(i1-1,j1-1)f1(i1-1,j1)f1(i1-1,j1+1)f1(i1-1,j1+2)f1(i1,j1-1)f1(i1,j1)f1(i1,j1+1)f1(i1,j1+2)f1(i1,j1-1)f1(i1+1,j1)f1(i1+1,j1+1)f1(i1+1,j1+2)f1(i1+2,j1-1)f1(i1+2,j1)f1(i1+2,j1+1)f1(i1+2,j1+2)]；g1(i,j)=(A*B*C)；
end end
g=uint8(g1)；%将矩阵转换成8位无符号整数 imshow(g)；
title（'自编双三次插值法放大8倍图像'）；
其中子函数sw代码如下： function A=sw(w1)w=abs(w1)；if w=0 A=1-2*w^2+w^3;elseif w>=1&&w<2 A=4-8*w+5*w^2-w^3;else
A=0;end
与MATLAB自带函数相比，以上手工编写的MATLAB代码只能完成黑白图像输出，且运行时间远比MATLAB自带函数的运行时间长。手工编写双三次插值算法MATLAB代码的运行结果如图6所示：
图6
2，其他算法简介：
传统的图像放大方法有重复放大线性放大和高次多项式插值放大。重复放大最简单，但会产生明显的方块效应线性放大消除了方块效应，但会造成图像的模糊 高次多项式插值放大效果较好，但运算复杂。由于传统方法的固有缺陷，诞生了新一代图像放大方法，主要有小波放大、邻域交换内插和分形放大等。
下面简单介绍一下增强系数小波放大算法： 算法示意图如图7所示：
图7 通过二维离散小波变换，经分析高通滤波器和分析低通滤波器，可将一幅分辨率为p的二维图像分解为分辨率为p/2的离散逼近信号A1和水平、垂直、对角三个细节信号H1、V1、D1。这四个分量都只有原图像大小的1/4。之后又可以对A1进行同样的分解如图7所示。这个过程可以一直重复下去。通过二维离散小波反变换，用相应的综合高通滤波器和综合低通滤波器可将各分量重构为原图像。
对于一个图像，低频成分包含了基本特征，即原图像的近似，高频成分反应其细节。基于此，我们将原图像作为低频成分A1，其他3个细节部分置0，进行小波重构，便可得到放大4倍的图像。但是由于能量守恒，放大后的图像能量分散会显得较暗。可以将原图像灰度值矩阵乘2，再进行上述变换，便可解决这一问题。小波分解重构是一种全局运算，不会造成重复放大中的方块效应，同时较好地保持图像边缘的清晰。
%function [] = UseEMD(x,t)
function [] = UseEMD()
%UNTITLED2 Summary of this function goes here %Detailed explanation goes here
N = 39;% # of data samples
t = linspace(0,1,N)；
x=[20.34
21.25
20.62
19.3
17.6
15.68
15.46
16.27
17.94
18.97
18.71
18.47
19.11
19.31
18.67
17.32
16.21
16.09
16.04
17.54
19.3
20.11
20.42
21.33
21.94
22.17
22.97
23.09
23.39
23.97
26.71
29.31
29.82
30.35
29.78
27.94
27.17
27.83
27.36
]；
y=fft(x)；
plot(y,x)；
[imf,ort,nbits] = emd(x)；
emd_visu(x,t,imf,1)；
%¾ùÖµµÄƽ·½
imfp2=mean(imf,2)。^2
%ƽ·½µÄ¾ùÖµ
imf2p=mean(imf.^2,2)
%¸÷¸öIMFµÄ·½²î
mse=imf2p-imfp2
%·½²î°Ù·Ö±È£¬Ò²¾ÍÊÇ·½²î¹±Ï×ÂÊ
mseb=mse/sum(mse)*100
%ÏÔʾ¸÷¸öIMFµÄ·½²îºÍ¹±Ï×ÂÊ
[mse mseb]
%¼ÆËãimfµÄÐÐÁÐάÊý
[m,n] = size(imf)
hx = hilbert(x)
xf = real(hx)；
xi = imag(hx)；
%¼ÆËã˲ʱÕñ·ù
A=sqrt(xf.^2 + xi.^2)；
figure(4)
plot(t,A)；title（'˲ʱÕñ·ù'）
%¼ÆËã˲ʱÏàλ
p=angle(hx)；
figure(5)；plot(t,p)；title（'˲ʱÏàλ'）
%¼ÆËã˲ʱƵÂÊ
%dt = diff(t)；
%dx=diff(p)；
%sp = dx./dt;
%figure(6)；
%plot(t(1:N-1)，sp)；title（'˲ʱƵÂÊ'）
%imf1µÄhilbert±ä»»
xn1 = hilbert（imf(1,:）)；
xr1 = real(xn1)；
xi1 = imag(xn1)；
%imf1µÄ˲ʱÕñ·ù
A1=sqrt(xr1.^2+xi1.^2)；
figure(7)；
subplot(2,1,1)；
plot(t,A1)；
xlabel（'ʱ¼äª¡tª¢'）
；ylabel（'˲ʱÕñ·ù'）；
title（'imf1'）%imf2µÄHilbert±ä»»
xn2=hilbert（imf(2,:）)；
xr2=real(xn2)；
xi2=imag(xn2)；%imf2µÄ˲ʱÕñ·ù
A2=sqrt(xr2.^2+xi2.^2)；
subplot(2,1,2)；
plot(t,A2)；
xlabel（'ʱ¼äª¡tª¢'）；
ylabel（'˲ʱÕñ·ù'）；
title（'imf2'）%imf1µÄ˲ʱÏàλ
P1=atan2(xi1,xr1)；
figure(8)；
subplot(2,1,1)；
plot(t,P1)；
xlabel（'ʱ¼äª¡tª¢'）；
ylabel（'˲ʱÏàλ'）；
title（'imf1'）%imf2µÄ˲ʱÏàλ
P2=atan2(xi2,xr2)；
subplot(2,1,2)；
plot(t,P2)；
xlabel（'ʱ¼äª¡tª¢'）；
ylabel（'˲ʱÏàλ'）；
title（'imf2'）%imf1˲ʱƵÂÊ
xh1=unwrap(P1)；
fs=1000;
xhd1=fs*diff(xh1)/(2*pi)；
figure(9)；
subplot(2,1,1)；
plot(t(1:38)，xhd1)；
xlabel（'ʱ¼äª¡tª¢'）；
ylabel（'˲ʱƵÂÊ'）；
title（'imf1'）%imf2˲ʱƵÂÊ
xh2=unwrap(P2)；
fs=1000;
xhd2=fs*diff(xh2)/(2*pi)；
subplot(2,1,2)；
plot(t(1:38)，xhd2)；
xlabel（'ʱ¼äª¡tª¢'）；
ylabel（'˲ʱƵÂÊ'）；
title（'imf2'）%¼ÆËãHHTµÄʱƵÆ×
[A, fa, tt] = hhspectrum(imf)；
if size(imf,1)>1
[A,fa,tt] = hhspectrum（imf(《www．》1:end-1, :）)；
else [A,fa,tt] = hhspectrum(imf)；
end
[E,tt1]=toimage（A,fa,tt,length(tt)）；
disp_hhs(E,[]，fs)%¶þάͼÐÎÏÔʾHHTÊÓƵÆת£EÊÇÇóµÃµÄHHTÆ×for i=1:m-1
faa=fa（i,:）；
[FA,TT1]=meshgrid(faa,tt1)；%ÈýάͼÏÔʾHHTʱƵͼfigure(11)；
surf(FA,TT1,E)
title（'HHTʱƵÆ×ÈýάÏÔʾ'）
end
%»-±ß¼ÊÆ×
E=flipud(E)；
for k=1:size(E,1)
bjp（k)=sum(E(k,:）)*1/fs;
end
f =(0:N-3)/N*(fs/2)；
figure(12)
plot(f,bjp)；
end
注意：matlab中需加载instfreq.m文件，从网上可下到
最近矩阵分析老师出了一道题目作为作业，是一道程序题，题目是对A矩阵做LU分解，要求能对A实现PA=LU的分解，并最终输出L，U，P矩阵。
先来解读下题目，寥寥几句话，里面囊括的信息量却不少，然后这些都得自己去琢磨。首先对A矩阵能做LU分解，即能把A分解成这种形式A=L*U（U是上三角矩阵，是由A矩阵经过高斯消元后得到的，L是下三角矩阵，其对角线全为1，其他非零元素为在消去（i,j）位置元素过程中主元所乘的系数），条件有3，一是矩阵A必须为方阵，A如果不是方阵，就不要想着对它做LU分解啦，这是基本条件，牢记啊！二是矩阵A必须可逆，换种说法就是A必须为非奇异矩阵，这两种说法是等价的，而这又等价于A是满秩的，A是满秩又等价于A的行列式值非0，好绕，矩阵就是这样，很多定理其实是等价的，但是你得记住，不然在推导一些定理或公式的时候会犯一些基本的常识性错误。三是矩阵A在高斯消元过程中必须没有出现0主元，也就是说只有在对A进行高斯消元过程中没有出现进行交换这种情况下，A才能分解成L*U这种形式，如果对A进行高斯消元，中间某一步出现0主元，需要进行行交换了，这种情况下就不要想对A进行LU分解啦，因为不满足条件3啊！那么问题来了，假如出现了有0主元这种情况，我又想对A进行LU分解，那应该怎么办？
这就引出了带行变换的LU分解，也就是本文的主题。根据书上的定理，对任意一个非奇异矩阵（等价于可逆矩阵）都存在一个置换矩阵P使得P*A可以分解成L*U这种形式，即PA=LU。想想其实这定理也是不言自明的，刚才A不是要进行行变换才能继续高斯消元吗？而LU分解前提又是高斯消元过程中不能出现行交换，那好，我事先对A矩阵在高斯消元过程中需要交换的行给交换掉，形成一个新的矩阵B，那我对B高斯消元那肯定就不会出现需要行交换的情况，这就满足了LU分解第三个条件了，这样B不就可以进行LU分解了吗？是的，PA=LU这种形式的LU分解采用的就是这种思想。那么现在的问题是，怎么在代码中实现对A矩阵的LU分解，并输出P，L，U矩阵呢？我在网上搜了一下，发现结果大都不尽如人意，大多数程序吧只能说做A=LU这种形式的分解，一旦说A不满足条件3，那就死翘翘了，这种程序先不论其能否运行成功，结果是否正确，其鲁棒性也太差了！用个时髦点的词来说就是太low了！通用性太差了！不光如此，代码也没什么注释，可读性很差，让人怀疑是不是写给别人看的，尤其像我这种编程渣渣，看个代码费老半天劲都看不懂说什么。于是，我决定按自己的想法来走。
首先从最简单的情况考虑，这也是我们做研究、做学术、做工程必须要时刻牢记心中的一点，很多人喜欢一上来就把所有问题、把最复杂的情况、把方方面面都给考虑到，然后再开始实现他的想法，我自己也有这个习惯，但是，这并不是一个好习惯，一上来就好高骛远、就想着高大上这本质上是一种急功近利的表现，那样的话你会陷入到各种各样的技术细节当中，你会想半天却仍然写不出半点实质性的东西出来，所以最好的办法是，先考虑最简单、最核心的情况，这样不仅大大降低问题的复杂度，同时也为将来进一步扩展程序、解决更复杂的情况打下了一个坚实的基础。
在这个例子中，最简单的情况就是矩阵A在高斯消元过程中不需要进行行交换，也就是说A可以分解成A=L*U这种形式。这种情况下，代码如下。
function LUDecomposition(A,n)%A is the square matrix,n is the order of A
L=eye(n)；%Let the L matrix be an identity matrix at first for i=1:n-1
for j=i+1:n
L(j,i)=A(j,i)/A(i,i)；
A（j,:）=A（j,:）-（A(j,i)/A(i,i)）*A（i,:）；end end
U=A %A becomes U matrix after Gauss elimination L
可以试着令A=[2 2 2;4 7 7;6 18 22]，调用函数获得L矩阵为[1 0 0;2 1 0;3 4 1]，U矩阵为[2 2 2;0 3 3;0 4 4]，用笔验算下，这个结果是正确的。代码运行结果如图所示：
这部分代码的主要思想是这样的，矩阵A的阶次为n的话，A在高斯消元后有n个非零主元。在消元过程中，A共需要消掉n-1个主元下面所有的元素，注意，第n个主元已经是矩阵的最后一个元素了，它的下面和右边都没有其他元素了，所以不存在说对第n个主元下面所有元素消去的情况。这就获得了我们代码的第一个for循环，从第1行主元开始消元，一直到第n-1行主元。而在获得每一行主元过程中，需要对该行主元下面所有元素都消去，假如现在要获得第i行主元的话，就是说要对该主元所在列的第i+1行到第n行元素都消掉，那么这就获得了我们代码的第二个for循环，从消去第i+1个元素开始一直到第n个元素。前文说过，消掉第（j，i）个位置元素过程中，主元所乘系数就是L矩阵第（j，i）位置的元素，所以有L（j,i)=A(j,i)/A(i,i)。然后的话，就是把A矩阵第j行减去第j行乘以L（j，i），这样就可以消掉第（j，i）个元素了，就是这行代码A(j,:）=A（j,:）-（A(j,i)/A(i,i)）*A（i,:）。最后，执行完两层for循环后，A矩阵就成为了U矩阵，L矩阵也从最初的单位阵成了L矩阵。
好了，我们已经实现了最简单的情况了，下面考虑复杂点的情况，就是说对A进行PA=LU这种形式的分解。假如A在消元过程中出现0主元了，那么怎么办？很简单，只需要从该0主元下面所有元素中找到一个非0元素，然后将其所在的行与该0主元所在的行进行交换就行了，不要忘了，对A矩阵两行进行了交换，对应到P矩阵中的操作是相应的两行也要进行交换，因为我们是通过P矩阵两行交换后然后左乘A矩阵使得A矩阵两行进行交换的。A矩阵交换第i行和第k行元素对应到L矩阵中相应两行的消元系数也应该交换位置，就是说L矩阵的第i行和第k行元素也要交换位置，当然，主对角线上的1是不需要交换的，因为他们并不是消元系数。交换完成后，继续执行消元操作，其步骤和上面考虑的最简单的情况就是一样的了。就这样，我们就实现了PA=LU这种形式的分解。令A=[1 2-3 4;4 8 12-8;2 3 2 1;-3-1 1-4]，代入函数运算得L矩阵为[1 0 0 0;-3 1 0 0;2-0.2 1 0;2-0.2 1 0;4 0 3.75 1]，U矩阵为[1 2-3 4;0 5-8 8;0 0 6.4-5.4;0 0 0-3.75]，P矩阵为[1 0 0 0;0 0 0 1;0 0 1 0;0 1 0 0]，用笔验算下，结果与函数运行结果是一致的，当然了，这个函数我只是代了3，4个不同的A矩阵进去而已，可能样本数量不够多，但目前来说我觉得应该没什么问题了，如果有问题欢迎反馈给我。这部分代码如下：
function AdvanceLUDecomposition(A,n)%A is the square matrix,n is the order of A,A must be invertible
D=A;%Store matrix A in D,for later use
L=zeros(n)；%Let the L matrix be an zero matrix at first
P=eye(n)；%Let the permutation matrix be a identity matrix at first for i=1:n-1 for j=i+1:n
if A(i,i)==0 %A zero pivot appears on(i,i)position,we need to find a nonzero entry below it to be the new pivot,with row exchange for k=n:-1:i+1 %find a nonzero entry below the(i,i)entry in the i column,start from the last row
if A(k,i)~=0 %We have found a nonzero entry,to choose it as the new pivot,we need row exchange ki
L（[i k]，:）=L（[k i]，:）；%Permute i and k row in L matrix A（[i k]，:）=A（[k i]，:）；%Permute i and k row in A matrix P（[i k]，:）=P（[k i]，:）；%Permute i and k row in P matrix break;end end end
L(j,i)=A(j,i)/A(i,i)；
A（j,:）=A（j,:）-（A(j,i)/A(i,i)）*A（i,:）；end end
U=A %A becomes U matrix after Gauss elimination
L=L+eye(n)%All entries on the diagonal of L matrix must be 1 P %output the permutation matrix
B=L*U %verify if the product of L and U equals to P*A
C=P*D %D is the original A matrix,check it out in row 2 %If B equals C,then it means the algorithm works correctly
%some key points and theroms about LU factorization
%Theorem 1 A nonsigular matrix Anxn possesses an LU factorization if and %only if a nonzero pivot does not emerge during row reduction to upper %triangular form with type III operations.%Theorem 2 For each nonsigular matrix A,there exist a permutation matrix P
%such that PA possesses an LU factorization PA=LU
%Remember,the concept of nonsingular matrix is for square matrix,it means %that the determinant is nonzero,and this is equivalent that the matrix has %full-rank
%Based on these conditions,the first thing about the matrix A on which we
%conduct LU factorization is that A must be a square matrix.The second %thing is A must be invertible,which is equal to the statement that A is %non-singular
代码运行结果如图所示：
最后补充一点，为什么要进行LU分解呢？这个问题很关键，很多人也许并不关注这个问题，我们学习很多时候都是只关注实现方法，却并不关心它存在的意义，这种学习是永远无法深入的，只能是停留在表面上，学习就应该多问为什么，多质疑这个东西存在的价值，存在的意义有多大，这样才能促使你去深入了解这个方法的优点和缺点，从而改进、完善它。简单点来说就是LU分解大大降低了算法复杂度，我们求解一个方程组Ax=b的时候，一般来说无非就两种方法，要么是高斯消元法，要么是先求A的逆矩阵，然后再乘以b获得x，而第二种方法比第一种方法要复杂并且限制更多，所以一般是用高斯消元法。高斯消元法解一个方程组算法复杂度是(n^3)/3，并且每获得一个新的b，要接x，都得执行复杂度为(n^3)/3的操作。而LU分解有什么好处呢？在第一次LU分解的时候，也就是说获得L和U的时候，其算法复杂度其实也是n^3,但是，一旦我们获得了L和U矩阵后，每次我们获得一个新的b要求对应的解x，算法复杂度就会大大降低，粗略来说就是n^2，把复杂度降低了一个级别，对于大型系统来说，这是非常了不起的一个改进，运算性能会大大提升。而实际应用中，这样的方法也是非常有意义的，实际系统中，A矩阵相当于系统里的各种滤波和变换操作，x相当于系统的输入，b相当与系统的输出，我们一般是获得了输出b，然后想求得输入x，只要系统不变，那么知道b，又知道了L和U矩阵，我们只需要对每一个新的b执行n^2次乘法/除法和n^2-n次加法/减法就可以获得b对应的输入x了，这是多么了不起的一个性能改进！正因为这样，LU分解在实际应用中用的也是非常广泛。
写完这篇文章，不知道矩阵分析老师会不会看到呢？不知道他以后还会不会出这道题呢？假如还出这道题的话，希望我这篇文章能对还在苦苦寻找源代码的各位师弟师妹们能起到一点小小帮助，当然了，这也只是一个抛砖引玉的方法，希望各位看官能有更好的答案，请不吝赐教！
9加几
例：9
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=1108、7、6加几
例：8
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=
+
=
十几减9
=
=
=
=
=
=
=
=
十几减8、7、6
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
=
数学与信息科学学院
实 验 报 告
课程名称： 《 数值计算方法 》 实验名称： 数值积分 实验类型： 验证性■综合性□设计性□ 实验室名称： 数学实验室 班级学号： 090721 学生姓名：
任课教师（教师签名）： 成 绩：
实验日期： 2012-4-27
一、实验目的及题目
实验目的：
南昌航空大学数学与信息科学学院实验报告
熟悉线性方程组求解原理，运用列主消元高斯消去法，LU分解法及Jacobi迭代法与Gauss-Seidel迭代法丢出线性方程的解 实验题目：
1、用列主元消去法解方程组：
x1x23x442xxxx1123
4
3xxx3x32341x12x23x3x44
2、用LU分解法界方程组Axb，其中
48240124122412124
A，b
062022626216
3、分别用jacobi迭代法和gauss-seidel迭代法求解方程组：
10x1x22x3118xx3x1123
4
2xx10x6231x13x2x311x425
二、实验原理、程序框图、程序代码等
实验原理：
1、列主元高斯消去原理：
每次消去中，选取每列的绝对值最大的元素作为消去对象并作为主元素。然后换行使之变到主元位子上，在进行消元计算。设A(k)xb(k)，确定第k列主元所在位置ik，在交换ik行和k行后，在进行消元。根据矩阵理论，交换ik和k两个方程的位置，列主元素的消去过程相当于对交换后的新矩阵进行消元，即
LkIk,ikAkA(k1)
同时，右端向量b(k)变化为 LkIk,ikbkb(k1)
2、LU分解原理：
第 页
南昌航空大学数学与信息科学学院实验报告
直接三角分解：先将A分解为上三角矩阵U和下三角矩阵L（A=LU），则原为题等价于求解两个方程组：Lyb,Uxy
3、迭代法
（1）jacobi迭代法（又称简单迭代法）考虑n阶线性代数方程组
a11x1a12x2。.。.。a1nxnb1axax。.。.。axb2112222nn
2
an1x1an2x2。.。.。annxnbn其矩阵形式为
Axb
设该方程组的系数矩阵A非奇异且aii0(i1,2,。.。.。，n)，可将A分解为：
ADLU
其中Ddiag[a11,a22,。.。.。，ann]；
0a21La31an1 0a31an10an10a120，U=－0a13a230a1na2n
an1,n0然后化为如下等价形式
xD1(LU)xD1b
简记为
xBJxfJ
选取初始向量x(0)x1,x2(0)(0)xn(0)T，通过上式可得到线性代数方程组的迭代格式
x（k1)BJx(k)fJ,(k0,1,2）
其分量形式为
xi即（k1)n1k(biaijxj)，(i1,2,n)(k0,1,2）aiii1j1第 页
南昌航空大学数学与信息科学学院实验报告
(k1)1(k)(k)(k)x(baxaxax)11221331nn1a111(k1)(k)(k)(k)x(baxaxax)222112332nna 221(k1)(k)(k)(k)x(baxaxax)nnn1n2n,n1n112ann
以上计算过程称迭代法，矩阵BJ陈为jacobi迭代法的迭代矩阵。
（2）Gauss-Seidel迭代法 将jacobi迭代格式改为如下形式
xi（k1)i1n1(k)(k)(biaijxjaijxj)，(i1,2,n)(k0,1,2）aiii1ji1其矩阵形式为
x(k1)(DL)1Ux(k)(DL)1b
于是得Gauss-Seidel迭代公式为
x（k)BGx(k)fG,(k0,1,2）
程序代码： 第一题程序代码： function GEpiv（A,b)[m,n]=size(A)；nb=n+1;Ab=[A b]；for i=1:m-1 [pivot,p]=max(abs(Ab(i:n,j)）)；ip=p+i-1;if ip~=i Ab（[i ip]，:）=Ab（[ip i]，:）；end
第 页 称BG为Gauss-Seidel迭代法的迭代矩阵。
南昌航空大学数学与信息科学学院实验报告
pivot=Ab（i,i)；for k=i+1:m Ab(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb)；end end x=zeros(n,1)；x(n)=Ab(n,nb)/Ab(n,n)；for i=n-1:-1:1 x(i)(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n,1)）/Ab（i,i)；end for k=1:n fprintf（'x[%d]=%fn'，x(k））；end 运行结果为 输入
A=[1 1 0 3;2 1-1 1;3-1-1 3;-1 2 3-1]；b=[4;1;-3;4]；得结果为 x =
-1.3333 2.3333-0.3333 1.0000
第二题代码：
function x=lufj(A,b)%lu分解求值函数 [l u p]=lu(A)；y=lx(l,b)；x=us(u,y)；end
第 页
南昌航空大学数学与信息科学学院实验报告
function y=lx(A,b)%求系数矩阵为下三角方程函数 [n,m]=size(A)；y(1)=b(1)/A(1,1)；for i=2:n s=0;for k=1:i-1 s=s+A(i,k)*y(k)；end y(i)=(b(i)-s)/A(i,i)；end end
function x=us(A,b)%求系数矩阵为上三角方程函数 [n,m]=size(A)；x(n)=b(n)/A(n,n)；for i=0:n-2 s=0;for k=0:i s=s+A(n-i-1,n-k)*x(n-k)；end x(n-i-1)=(b(n-i-1)-s)/A(n-i-1,n-i-1)；end end 输入矩阵
A=[48-24 0 12;-24 4 12 12;0 6 20 2;-6 6 2 16]；b=[4;4;-2;-2]；结果为 x =
第 页
南昌航空大学数学与信息科学学院实验报告
0.5212 1.0055-0.3757-0.2597
第三题代码：（1）Jacobi迭代：
function x = agui_jacobi（a,b)n=length(b)；N=100;e=1e-4;x0=zeros(n,1)；x=x0;x0=x+2*e;k=0;d=diag(diag(a)）；l=-triu(a,-1)；u=-triu(a,1)；while norm(x0-x,inf)>e&k （2）Gauss-seidel迭代 function x=gau(A,b,x)[d u l]=fenjie(A)；B=inv(d-l)*u;f=inv(d-l)*b; 第 页 南昌航空大学数学与信息科学学院实验报告 for i=1:9 y=x;x=B*y+f;end function [d u l]=fenjie(A)[n m]=size(A)；for i=1:n for j=1:n d(i,j)=0;end end for i=1:n d(i,i)=A(i,i)；end for i=1:n for j=1:i l(i,j)=-A(i,j)；u(j,i)=-A(j,i)；end for k=i:n l(i,k)=0;u(k,i)=0;end end（1） 输入矩阵：a=[10-1 2 0;0 8-1 3;2-1 10 0;-1 3-1 11]；>>b=[-11;-11;6;25]；>>x = agui_jacobi(a,b)结果为 第 页 南昌航空大学数学与信息科学学院实验报告 x=-1.4674-2.3587 0.6576 2.8424 （2）A=[10-1 2 0;0 8-1 3;2-1 10 0;-1 3-1 11]；b=[-11;-11;6;25]；x=[-1.3;-2.1;0.55;2.66]；gau(A,b,x)结果为 x=-1.4674-2.3587 0.6576 2.8424 四、实验中存在的问题及解决方案 命令需多次修改，如Gauss迭代中，m文件中因end的不匹配使得程序无法运行处结果，则需细心调试，去除多余的end，才能使程序正常运行。 五、心得体会 在线性代数中，从理论上解决了线性方程组的求解问题，但常用的高斯消元法有时会导致结果产生严重的偏差，因此需要研究数值解法。Jacobi迭代与Gauss-Seidel有一个共同特点：新的近似解是已知近似解的线性函数，并且新的近似解只与已知近似解有关。 编写程序是一个繁杂的过程，需要不断调试，并且善于发现程序中的细小问题，才能够正常运行处程序。 第 页