matlab中复合中点式程序,《现代数值计算》Matlab程序整理（23页）-原创力文档

第2 章 线性方程组的直接解法

1.1 写一个求解线性代数方程x b 的列主元素高斯消去法程序，在程序入口输入n 的

1

值，其中ij ，bj ln j ，1 i , j  n

i  j 1

解：

function x=gauss(n)

for i=1:n

for j=1:n

A(i,j)=1/(i+j-1);

end

b(i,1)=log(i);

end

%Gauss Elimination

for i=1:n-1

[r,p]=max(abs(A(1:n,i)));

p=p+i-1;

A([p,i],:)=A([i,p],:);

for j=i+1:n

r -A(j,i)/A(i,i);

A(j,:)=A(j,:)+r*A(i,:);

b(j)= b(j)+r*b(i);

end

end

x=zeros(n,1);

for k=n:-1:1

x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);

end

1

1.2 写一个求解线性代数方程A=[1,13,-2,-34;2,6,-7,-10;-10,-1,5,9;-3,-5,0,15];b=[13,-22,14,-36]';

的列主元素高斯消去法程序

解：

function x=gaussel(A,b)

n=length(b);

for i=1:n-1

[r,p]=max(abs(A(i:n,i))); %r 为最大值，p 为所在行

p=p+i-1;

if p>i

t=A(i,:);A(i,:)=A(p,:);A(p,:)=t;

t=b(i);b(i)=b(p);b(p)=t;

end

for j=i+1:n

r -A(j,i)/A(i,i);

A(j,:)=A(j,:)+r*A(i,:);

b(j)= b(j)+r*b(i);

end

end

x=zeros(n,1);

x(n)=b(n)/A(n,n);

for k=(n-1):-1:1

x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);

end

>> A=[1,13,-2,-34;2,6,-7,-10;-10,-1,5,9;-3,-5,0,15];

>> b=[13,-22,14,-36]';

>> x=gaussel(A,b)

x

14.0669

29.9869

17.9950

10.4090

2

2. LU 分解

function [L,U]=myLU(A) %实现对矩阵A 的LU 分解，L 为下三角矩阵

[n,n]=size(A);

L=zeros(n,n);

U=zeros(n,n);

for i=1:n

L(i,i)=1;

end

for k=1:n

for j=k:n

U(k,j)=A(k,j)-sum(L(k,1:k-1)*U(1:k-1,j));

end

for i=k+1:n

L(i,k)=(A(i,k)-sum(L(i,1:k-1)*U(1:k-1,k)))/U(k,k);

end

end

3