数值线性代数实验

数值线性代数实验

题目：数值线性代数

专业：信息与计算科学班级：班姓名：

山东科技大学

2013年 1 月16日

实验报告说明

学院：信息学院专业：信息班级10-2 姓名：

一、主要参考资料：

（1）《Matlab数值计算-案例分析》北京航空出版（2）《Matlab数值分析》机械工业出版

二、课程设计应解决的主要问题：

（1）平方根

（2）QR方法

（3）最小二乘法

三、应用软件：

（1）Matlab7.0

（2）数学公式编辑器

四、发出日期：课程设计完成日期：

指导教师签字：系主任签字：

指导教师对课程设计的评语

指导教师签字：

年月日

一、问题描述

先用你所熟悉的计算机语言将平方根和改进的平方根法编成写通用的子程序，然后用你编写的程序求解对称正定方程组b x =A ，其中 （1）b 随机的选取，系数矩阵位100阶矩阵

⎥⎥

⎥

⎥⎥⎥

⎥⎥⎥⎦

⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡1011101110111011101110

（2）系数矩阵为40阶Hilbert 矩阵，即系数矩阵A 的第i 行第j 列元素为

11-+=j i a ij ，向量b 的第i 个分量为∑=-+=n

j i j i b 11

1

。

二、分析与程序

1. 平方根法函数程序如下：

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)

求解Ax=b

L=zeros(n,n);

D=diag(n,0);

S=L*D;

for i=1:n %L的主对角元素均为1

L(i,i)=1;

end

for i=1:n

for j=1:n

if (eig(A)<=0)|(A(i,j)~=A(j,i))

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));

end

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

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 %LDL'分解

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

3.调用函数解题：

clear;clc;

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

for i=1:n

for j=1:n