辽宁工程技术大学上机实验报告

 

第二篇：数值代数上机实验报告

试验项目名称：平方根法与改进平方根法

实验内容：先用你熟悉的计算机语言将平方根法和改进平方根法编写成通用的子程序，然后用你编写的程

序求解对称正定方程组

Ax=b，其中，

A=[10

1 10 1

…

1 10 1

1 10]100*100

b随机生成，比较计算结果，评论方法优劣。

实验要求：平方根法与改进的平方根的解法步骤；存储单元，变量名称说明；系数矩阵

与右端项的生成；结果分析。

实验报告 姓名：罗胜利 班级：信息与计算科学0802 学号：实验一、平方根法与改进平方根法

先用你所熟悉的计算机语言将平方根法和改进的平方根法编成通用的子程序，然后用你编写的程序求解对称正定方程组AX=b，其中 系数矩阵为40阶Hilbert矩阵，即系数矩阵A的第i行第j列元素为=，

向量b的第i个分量为=.

平方根法函数程序如下：

function [x,b]=pingfanggenfa(A,b)

n=size(A);

n=n(1);

x=A^-1*b; %矩阵求解 disp('Matlab自带解即为x');

for k=1:n

A(k,k)=sqrt(A(k,k));

A(k+1:n,k)=A(k+1:n,k)/A(k,k);

for j=k+1:n;

A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k);

end

end

for j=1:n-1

b(j)=b(j)/A(j,j);

b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);

end

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

A=A';

for j=n:-1:2

b(j)=b(j)/A(j,j);

b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);

end

b(1)=b(1)/A(1,1);

disp('平方根法的解即为b');

end

function [x]=ave(A,b,n)

L=zeros(n,n);

D=diag(n,0);

S=L*D;

for i=1:n

L(i,i)=1; %Cholesky分解 %前代法 %回代法 %用改进平方根法求解Ax=b %L为n*n矩阵 %D为n*n的主对角矩阵 %L的主对角元素均为1

end

for i=1:n

for j=1:n %验证A是否为对称正定矩阵

if (eig(A)<=0)|(A(i,j)~=A(j,i)) %A的特征值小于0或A非对称时，输出wrong disp('wrong');break;end

end

end

D(1,1)=A(1,1);

for i=2:n

for j=1:i-1

S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');

L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);

end

D(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');

end

y=zeros(n,1);

x=zeros(n,1);

for i=1:n

y(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i);

end

for i=n:-1:1

x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n)); %将A分解使得A=LDLT % x，y为n*1阶矩阵 %通过LDy=b解得y的值 %通过LTx=y解得x的值

end

改进平方根法函数程序如下：

function b=gaijinpinfanggenfa(A,b)

n=size(A);

n=n(1);

v=zeros(n,1);

for j=1:n

for i=1:j-1

v(i)=A(j,i)*A(i,i);

end

A(j,j)=A(j,j)-A(j,1:j-1)*v(1:j-1);

A(j+1:n,j)=(A(j+1:n,j)-A(j+1:n,1:j-1)*v(1:j-1))/A(j,j); end B=diag(A);

D=zeros(n);

for i=1:n

D(i,i)=B(i);

A(i,i)=1;

End

A=tril(A);

for j=1:n-1

b(j)=b(j)/A(j,j);

b(j+1:n)=b(j+1:n)-b(j)*A(j+1:n,j);

end

b(n)=b(n)/A(n,n); A=D*(A');

for j=n:-1:2

b(j)=b(j)/A(j,j);

b(1:j-1)=b(1:j-1)-b(j)*A(1:j-1,j);

end

b(1)=b(1)/A(1,1); disp('改进平方根法解得的解即为b');

end

调用函数解题：

clear;clc;

n=input('请输入矩阵维数：');b=zeros(n,1); A=zeros(n);

for i=1:n

for j=1:n

A(i,j)=1/(i+j-1); %LDL'分解 %得到L和D %前代法 %回代法

b(i)=b(i)+1/(i+j-1);

end

end %生成hilbert矩阵

[x,b]=pingfanggenfa(A,b)

b=gaijinpinfanggenfa(A,b)

运行结果：

请输入矩阵维数：40

Warning: Matrix is close to singular or badly scaled.

Results may be inaccurate. RCOND = 6.570692e-020. > In pingfanggenfa at 4

In qiujie at 10

Matlab自带解即为x

平方根法的解即为b

x =

1.6035

8.9685

0.8562

1.0195

0.9375

-50.2500

-3.0000

-16.0000

24.0000

-49.5000

-30.0000

39.0000

22.0000

-64.0000

-12.0000

2.0000

10.2500

-10.5000

-1.0000

-10.8750

83.0000

46.0000

-98.0000

12.0000

-69.0000

68.0000

21.0000

17.0000 -50.7188 -8.7500 -8.0000 112.0000 6.0000 -68.7500 22.0000 44.0000 -28.0000 8.0000 -44.0000 12.0000

b =

1.0e+007 *

0.0000 -0.0000 0.0001 -0.0004 -0.0014 0.0424 -0.2980

1.1419 -2.7335

4.2539 -4.3018

2.7733 -1.1989 0.5406 -0.3688 0.3285 -0.4438 0.4621 -0.2513 0.0565 0.0000 -0.0051 0.0071 -0.0027 -0.0031

0.0036

-0.0019

0.0009

0.0002

-0.0002

-0.0006

0.0004

0.0001

-0.0002

0.0001

0.0000

-0.0000

0.0000

-0.0000

-0.0000

改进平方根法解得的解即为b

b =

1.0e+024 *

0.0000

-0.0000

0.0001

-0.0012

0.0139

-0.0954

0.4208

-1.2101

2.0624

-1.0394

-3.3343

6.2567

-0.2463

-7.4594

2.8030

3.6990

0.7277

-1.7484

-0.4854

-3.6010

0.2532

5.1862

-2.1299

1.4410

0.8738

-4.5654

1.0422

4.0920

-2.7764

-2.2148

-0.8953

0.3665

4.8967

1.0416

0.1281

-4.3387

-1.1902

-2.8334

8.4610

-3.6008

实验二、利用QR分解解线性方程组：

利用QR分解解线性方程组Ax=b，其中 A=[16 4 8 4;

4 10 8 4;

8 8 12 10;

4 4 10 12];

b=[32 26 38 30];

求解程序如下：

定义house函数：

function [v,B]=house(x)

n=length(x);

y=norm(x,inf);

x=x/y;

Q=x(2:n)'*x(2:n);

v(1)=1;v(2:n)=x(2:n);

if n==1

Q=0;

B=0;

else

a=sqrt(x(1)^2+Q);

if x(1)<=0

v(1)=x(1)-a;

else

v(1)=-Q/(x(1)+a);

end

B=2*v(1)^2/(Q+v(1)^2);

v=v/v(1);

end

end

进行QR分解：

clear;clc;

A=[16 4 8 4;4 10 8 4;8 8 12 10;4 4 10 12];

b=[32 26 38 30];

b=b';

x=size(A);

m=x(1);

n=x(2);

d=zeros(n,1);

for j=1:n

[v,B]=house(A(j:m,j));

A(j:m,j:n)=(eye(m-j+1)-B*(v')*v)*A(j:m,j:n); d(j)=B;

if j<m

A(j+1:m,j)=v(2:m-j+1);

end

end %QR分解 R=triu(A); %得到R D=A;

I=eye(m,n);

Q=I;

for i=1:n

D(i,i)=1;

end

H=tril(D);

M=H';

for i=1:n

N=I-d(i)*H(1:m,i)*M(i,1:m);

Q=Q*N;

end %得到Q

b=(Q')*b; %Q是正交阵 for j=n:-1:2

b(j)=b(j)/R(j,j);

b(1:j-1)=b(1:j-1)-b(j)*R(1:j-1,j);

end

b(1)=b(1)/R(1,1); %回带法

运行结果如下：

R =

18.7617 9.8072 15.7769 11.0864

0 9.9909 9.3358 7.5341

0 0 5.9945 9.8013

0 0 0 -0.5126

Q =

0.8528 -0.4368 -0.2297 -0.1709

0.2132 0.7916 -0.4594 -0.3417

0.4264 0.3822 0.2844 0.7689

0.2132 0.1911 0.8095 -0.5126

b=

1.00000000000000

1.00000000000001

0.999999999999988

1.00000000000001

实验三、Newton下山法解非线性方程组：

3x-cos（yz）-=0，

-81+sinz+1.06=0，

exp（-xy）+20z+=0；

要求满足数值解=满足或.

定义所求方程组的函数：Newtonfun.m

function F = Newtonfun(X)

F(1,1)=3*X(1)-cos(X(2)*X(3))-1/2;

F(2,1)=X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;

F(3,1)=exp(-X(1)*X(2))+20*X(3)+(10*pi-3)/3;

End

向量求导：Xiangliangqiudao.m

function J=xiangliangqiudao()

syms x y z

X=[x,y,z];

F=[3*X(1)-cos(X(2)*X(3))-1/2;X(1)^2-81*(X(2)+0.1)^2+sin(X(3))+1.06;exp(-X(1)*X(2))+20*X(3)+(10*pi-3)/3];

J=jacobian(F,[x y z]);

End

代值函数：Jacobi.m

function F=Jacobi(x)

F=[ 3,x(3)*sin(x(2)*x(3)), x(2)*sin(x(2)*x(3));2*x(1), -162*x(2)-81/5,cos(x(3));-x(2)/exp(x(1)*x(2)),-x(1)/exp(x(1)*x(2)),20];

End

方程组求解：

format long; %数据表示为双精度型 X1=[0,0,0]';

eps=10^(-8);

k=1;i=1;

X2=X1-Jacobi(X1)^(-1)*Newtonfun(X1);

while (norm(subs(X2-X1,pi,3.1415926),2)>=eps)&&(norm(Newtonfun(X1),2)>=eps)

if norm(Newtonfun(X2),2)<norm(Newtonfun(X1),2) %判断先后两次迭代的大小 X1=X2;

B=inv(Jacobi(X2));

C=Newtonfun(X2);

X2=X2-B*C;

i=i+1;

else

v=1/(2^k); %引入下山因子

X1=X2;

B=inv(Jacobi(X2));

C=Newtonfun(X2);

X2=X2-v*B*C;

k=k+1;

end

end

j=i+k-1 %迭代次数

X=X2 %输出结果

运行结果如下： j =

5

X =

0.500000000000000 -0.000000000000000 -0.523598775598299