matlab上机作业报告(计算初等反射阵,用Householder变换法对矩阵A

一、给定向量x≠0，计算初等反射阵Hk。 1．程序功能：

给定向量x≠0，计算初等反射阵Hk。 2．基本原理：

若x??x,x,?,x??R的分量不全为零，则由

???sign(x1)x2??u?x??e1?(x1??,x2,?,xn)?1?2 ??u??(??x1)?22?T?uu?1T?I??uu?H?I?22u2??确定的镜面反射阵H使得Hx???e?y；当(1?k?n)时，由

n?21/2??sign(x)((x))?kki?i?k??u(k)?(0,?,0,xk??k,xk?1,?,xn)T?Rn ????1(u(k))Tu(k)??(??x)=?(u(k))Tkkkk?k2??1(k)(k)T?Hk?I??ku(u)Tn有Hkx?(x1,x2,?,xk?1,??k,0,?,0)?R

算法：

(1)输入x，若x为零向量，则报错 (2)将x规范化，

M?max?x,x,?,x?

如果M＝0，则报错同时转出停机 否则xi?xiM,i?1,2,?,n (3)计算??x2，如果x1?0，则????

(4)???(??x1) (5)计算u?x,u(1)?x1?? (6)H?I??uu (7)y?(?M?,0,?,0)

?1T

(8)按要求输出,结束

3．变量说明：

x － 输入的n维向量；

n － n维向量x的维数；

M － M是向量x的无穷范数，即x中绝对值最大的一项的绝对值； p － Householder初等变换阵的系数ρ； u － Householder初等变换阵的向量U s x n

－ 向量x的二范数； － 输入的n维向量； － n维向量x的维数；

p － Householder初等变换阵的系数ρ； u － Householder初等变换阵的向量U

k － 数k，H*x=y，使得y的第k+1项到最后项全为零； 4．程序代码：

(1)

function [p,u]=holder2(x)

%HOLDER2 给定向量x≠0,计算Householder初等变换阵的p,u

%程序功能：函数holder2给定向量x≠0,计算Householder初等变换阵的p,u； %输入：n维向量x；

%输出：[p,u]。p是Householder初等变换阵的系数ρ， % u是Householder初等变换阵的向量U。 n=length(x); % 得到n维向量x的维数； p=1;u=0; % 初始化p,u； M=max(abs(x)); if M==0

% 得到向量x的无穷范数，即x中绝对值最大的一项的绝对值；

% 如果x=0，提示出错,程序终止;

disp('Error: M=0');

return; else

x=x/M; % 规范化 end;

s=norm(x); if x(1)<0 s=-s; end u=x;

% 求x的二范数

% 首项为负，s值要变号

% 除首项外，其余各项x,u相同 % 计算u的首项 % 计算p

u(1)=s+x(1); p=s*u(1); (2)

if n==1 u=0; end % 若x是1×1维向量，则u=0

function H=holderk(x,k)

%HOLDERK 给定向量x≠0,数k，计算初等反射阵Hk，使HkX=Y，其中Y的第k+1项到

最后项全为零；

%程序功能：函数holderk给定向量x≠0，数k，计算初等反射阵Hk，使HkX=Y，%程序功能：函数holder2给定向量x≠0,计算Householder初等变换阵的p,u； %输入：n维向量x,数k；

%输出：H。H是Householder初等变换阵，H*x=y，使得y的第k+1项到最后项全为零； %引用函数：holder2；

n=length(x); % 得到n维向量x的维数； if k>n %如果k值溢出，报错； disp('Error: k>n');

end

H=eye(n); % 初始化H，并使H(1:k,1:k)=I；

[p,u]=holder2(x(k:n)); % 得到计算Householde初等变换阵的系数ρ、向量U； H(k:n,k:n)=eye(n-k+1)-p\%u*u'; % 计算H(k:n,k:n)=I-p\%u*u'； 5．使用示例: 情形1: X为零向量 >> x=[0,0,0,0]'; >> H=holderk(x,1) Error: M=0 H =

1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 情形2： K值溢出： >> x=[1,2,3,4]';

>> H=holderk(x,5) Error: k>n 情形3： K值为1： >> x=[2,3,4,5]'; >> H=holderk(x,1) H =

-0.2722 -0.4082 -0.5443 -0.6804 -0.4082 0.8690 -0.1747 -0.2184 -0.5443 -0.1747 0.7671 -0.2911 -0.6804 -0.2184 -0.2911 0.6361 检验: >> det(H)

ans =

-1.0000 >> H*x ans =

-7.3485 0.0000 0.0000 0.0000 情形4： (1)K值为3： >> x=[4,3,2,1]'; >> H=holderk(x,3) H =

1.0000 0 0 1.0000 0 0 0 0 检验:

>> det(H)

ans = -1

>> H*x

ans =

4.0000 3.0000 -2.2361 0 (2)K值为2: >> x=[4,3,2,1]'; >> H=holderk(x,2)

0 0 0 0 -0.8944 -0.4472 -0.4472 0.8944

H =

1.0000 0 0 0 0 -0.8018 -0.5345 -0.2673 0 -0.5345 0.8414 -0.0793 0 -0.2673 -0.0793 0.9604

>> det(H) ans =

-1.0000 >> H*x

ans =

4.0000 -3.7417 0.0000

0

二、设A为n阶矩阵，编写用Householder变换法对矩阵A作正交分解的程序。 1．程序功能：

给定n阶矩阵A，通过本程序用Householder变换法对矩阵A作正交分解，得出A＝QR 2．基本原理：

任一实列满秩的m×n矩阵A，可以分解成两个矩阵的乘积，即A＝QR，其中Q是具有法正交列向量的m×n矩阵，R是非奇异的n阶上三角阵。 算法：

(1)输入n阶矩阵A

(2)对A?k:n,k?，求Househoulder初等反射阵的?k,u(3)计算上三角阵R，仍然存储在A

HkA(k)(k)。k?1,2,?,n?1

?A(k?1)

j?k,k?1,?,n

1?n?k?(k)?tj???uiaij??k?i?k? aij(k?1)?aij(k)?tjui?k?(4)计算正交阵Q

QQ(1)(k)?IHk?Q(k?1)

i?1,2,?,n

ti?1n?k?qj?k(k)(k)iju?k?j j?k,k?1,?,nqij(k?1)

?k?j?qij?tiu(5)按要求输出,结束

3．变量说明：

A － 输入的n阶矩阵，同时用于存储上三角阵R； n Q

－ 矩（方）阵A的阶数；

－ Q是具有法正交列向量的n阶矩阵；

p,u － 向量A（k:n,k），对应初等反射阵的ρ,u k，jj，ii － 循环变量；

t1 － 计算上三角阵R的系数tj； t2 － 计算正交矩阵Q的系数ti； 4．程序代码：

function [Q,A]=qrhh(A)

%QRHH 用Householder变换法对n阶矩阵A作正交分解A=QR；

%程序功能：函数qrhh用Householder变换法对矩阵A作正交分解A=QR； %输入：n阶矩阵A；

%输出：[Q,A]。Q是具有法正交列向量的n阶矩阵，

% A（即R）是非奇异的n阶上三角阵，仍用输入的矩阵A存储。 %引用函数：

% holder2；示例 [p,u]=holder2(x); [n,n]=size(A); %求矩（方）阵A的阶数； Q=eye(n);

%构造正交矩阵Q(1)=I；

for k=1:n-1

[p,u]=holder2(A(k:n,k)); %向量A（k:n,k），对应初等反射阵的ρ,u

for jj=k:n %计算上三角阵R（仍存贮于A） t1=dot(u,A(k:n,jj))/p; %利用向量内积求和 A(k:n,jj)=A(k:n,jj)-t1*u; end

for ii=1:n

%计算正交矩阵Q

t2=dot(u,Q(ii,k:n))/p; %利用向量内积求和 Q(ii,k:n)=Q(ii,k:n)-t2*u'; end end

5．使用示例: (1)A为3阶矩阵： >> A=[1 2 3; 2 3 0; 3 4 5] A =

1 2 3 2 3 0 3 4 5

>> [q,r]=qrhh(A) q =

-0.2673 0.8729 0.4082 -0.5345 0.2182 -0.8165 -0.8018 -0.4364 0.4082 r =

-3.7417 -5.3452 -4.8107 0 0.6547 0.4364 -0.0000 0.0000 3.2660 检验： >> q*r

ans =

1.0000 2.0000 3.0000 2.0000 3.0000 0.0000

3.0000 4.0000 5.0000

(2)A为4阶矩阵：

>> A=[1 2 3 4; 2 3 0 1; 3 4 5 6;1 6 8 0] A =

1 2 3 4 2 3 0 1 3 4 5 6 1 6 8 0

>> [q,r]=qrhh(A)

q =

-0.2582 0.0597 -0.2660 -0.9268 -0.5164 -0.1045 0.8434 -0.1049 -0.7746 -0.2688 -0.4662 0.3323 -0.2582 0.9556 -0.0222 0.1399 r =

-3.8730 0 0 0 检验： >> q*r

ans =

1.0000 2.0000 3.0000 1.0000

-6.7132 -6.7132 -6.1968 4.4647 6.4805 -1.4783 -0.0000 -3.3070 -3.0178 0.0000 0 -1.8187 2.0000 3.0000 4.0000 3.0000 -0.0000 1.0000 4.0000 5.0000 6.0000 6.0000 8.0000 0

数值求解正方形域上的Poisson方程边值问题

???2u?2u????2?2??f(x,y)?xy,0?x,y?1?y????x??u(0,y)?u(1,y)?u(x,0)?u(x,1)?1

用MATLAB语言编写求解此辺值问题的算法程序，采用下列三种方法，并比较三种方法的计算速度。1、用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子；2、用块 Gauss-Sediel迭代法求解线性方程组Au=f；3、（预条件）共轭斜量法。 用差分代替微分，对Poisson方程进行离散化，得到五点格式的线性方程组

4ui,j?ui?1,j?ui?1,j?ui,j?1?ui,j?1?hfij2?i,j?Nu1,j?uN?1,j?ui,1?ui,N?1?1，写成矩阵形式Au=f。其中 ?A212??IA???????IA33????I??v2????v3??u??????I????vAN?1,N?1??N??T2 1?i,j?N?1

?b2?b3f??????bN?4?1????14??Aii?????????1??????1??4?v2?(u2,2,u3,2,...,uN,2),v3?(u2,3,u3,3,...,uN,3),......,vN?(u2,N,u3,N,...,uN,N)22TTTTb2?h(f2,2,f3,2,...,fN,2)?(u2,1?u1,2,u3,1,...,uN,1?uN?1,2),b3?h(f2,3,f3,3,...,fN,3)?(u1,3,0,...,fN?1,3),......,bN?h(f2,N,f3,N,...,fN,N)?(u2,N?1?u1,N,u3,N?1,...,uN,N?1?uN?1,N)h?fi,j

2TTTT1N?xiyj,,取N?10,则h?0.1i,j?2,3,...,N

一 用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。

1. 基本原理：

Gauss-Seidel迭代法计算简单，但是在实际计算中，其迭代矩阵的谱半径常接近1，因此收敛很慢。为了克服这个缺点，引进一个加速因子（又称松弛因子）对Gauss-Seidel方法进行修正加速。

假设已经计算出第k步迭代的解（i=1，2，···，n），要求下一步迭代的解（i=1，2，···，n）。首先，用Gauss-Seidel迭代格式计算x(k?1)

? x(k?1)??a12x(k)?a13x(k)???a1nx(k)?b12n?1a11a113a11a11?

?(k?1)a23(k)a2n(k)a21(k?1)b2x??x?x???x??13n 2a22a22a22a22??? ?an1(k?1)an2(k?1)ann?1(k?1)bn?(k?1) x??x?x???x?12n?1?naaaannnnnnnn?

然后引入松弛因子，用松弛因子对和作一个线性组合。

xi(k?1)??xi(k?1)?(1??)xi(k)，i=1,2,…,n

将二者合并成为一个统一的计算公式：

(k?1)(k)(k)(k)?a11x1??(?a12x2???a1nxn?b1)?(1??)a11x1?(k?1)(k?1)(k)(k)(k) ??(?a21x1?a23x3???a2nxn?b2)?(1??)a22x2?a22x2? ??(k?1)(k)?ax(k?1)??(?ax(k?1)?ax(k?1)???axn?1?bn)?(1??)annxnnnnn11n22nn?1?

2. 算法

(1)Gauss-Seidel迭代法引入松弛因子w： (k?1)(k?1)(k)ui,j??ui,j?(1??)ui,j,五点格式即为：

du?ui,j?(k?1)

i?2,3,?,n?ui,j(k)?4(bi,j?ai?1,jui?1,j?ai,j?1ui,j?1?ai,j?1ui,j?1?ai?1,jui?1,j?ai,jui,j)(k?1)(k?1)(k)(k)(k)(2)计算步骤：

第一步：给松弛因子赋初值w=1.1～1.8,给场值u和场源b赋初值 第二步：用不同的w进行迭代计算。置error=0;

计算

ui,j(k?1)?ui,j?(k)?4(bi,j?ui?1,j?ui,j?1?ui,j?1?ui?1,j?4ui,j(k?1)(k?1)(k)(k)(k?1))

在计算机上采用动态计算形式

??(k?1)(k?1)(k)(k)(k?1) du?(bi,j?ui?1,j?ui,j?1?ui,j?1?ui?1,j?4ui,j)?4?

?ui,j?ui,j?du(i,j?2，3，???，n)?

如果|du|>error则error=|du|,如果error

第三步：比较不同的w的迭代次数，用kk存放最小迭代次数，用ww和uu存放相应

的w及u。

3． 程序

① [u,k]=SOR(u,b,w) %%%%%%%（被下面程序调用）

%输入场初值u0、场源b及松弛因子w，通过五点差分格式进行迭代运算，

%如果第k+1次的迭代结果与第k次的差小于精度，则可以近似认为第k+1次的迭代 %结果是精确解，然后返回迭代次数k和迭代解

function [u,k]=SOR(u,b,w) %输出迭代结果u，及迭代次数k m=length(u); %m为u的维数 h=1/(m-1); %h为步长 N=10000;e=0.0000001; %e为精度

for k=1:N %k为记录迭代次数 error=0; for j=2:m-1

for jj=2:m-1

sum=4*u(jj,j)-u(jj-1,j)-u(jj+1,j)-u(jj,j-1)-u(jj,j+1);

du=w*(h^2*b(jj,j)-sum)/4; %计算u的修正量 u(jj,j)= u(jj,j)+du; %修正u if error

end

if error<=e,break;end %判断是否达到精度 end

② [kk,ww,uu]=SOR_5dianchafen(n)

%用超松弛迭代法求解正方形域上的Poisson方程边值问题 %用5点差分格式求取泊松问题

%输入n,对x、y轴进行n等分；先确定场u的边界及场源b,在调用[u,k]=SOR(u,b,w)； %用不同w计算的迭代次数不同，用kk存放最小的迭代次数， %用ww和uu分别存放最佳松弛因子w和精确解

function [kk,ww,uu]=SOR_5dianchafen(n)

w=[1.1:0.1:1.8];m=length(w); %w为松弛因子 kk=1000; ww=0; %kk是最少迭代次数，ww是最松弛因子 h=1/n; %h步长

b=zeros(n+1,n+1); %计算场源b tic;

for i=2:n+1 for j=2:n+1

b(i,j)=(i-1)*(j-1)*h^2; end end

uu=zeros(n+1,n+1); u=zeros(n+1,n+1); %对u赋初值 u(1,1:n+1)=1;u(n+1,1:n+1)=1;u(1:n+1,1)=1;u(1:n+1,n+1)=1; mu=u; %初值mu以便不同的w计算 for i=1:m %用不同的w计算迭代

[u,k]=SOR(mu,b,w(i)); %调用[u,k]=SOR(u,b,w)，返回迭代次数及精确解 if kk>k, kk=k;ww=w(i);uu=u;end %把最少迭代次数付给kk,及其w,u赋给ww,uu end

t=toc %统计程序运算时间

4.计算结果： >> format short >> n=10;

>> [kk,ww,uu]=SOR_5dianchafen(n) t = 0.0310 kk = 48 ww = 1.6000 uu =

Columns 1 through 8

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

1.0000 1.0011 1.0022 1.0031 1.0039 1.0044 1.0047 1.0045

1.0000 1.0022 1.0042 1.0061 1.0076 1.0087 1.0091 1.0088

1.0000 1.0031 1.0061 1.0088 1.0110 1.0126 1.0133 1.0128

1.0000 1.0039 1.0076 1.0110 1.0138 1.0159 1.0168 1.0162

1.0000 1.0044 1.0087 1.0126 1.0159 1.0183 1.0194 1.0189

1.0000 1.0047 1.0091 1.0133 1.0168 1.0194 1.0208 1.0203

1.0000 1.0045 1.0088 1.0128 1.0162 1.0189 1.0203 1.0201

1.0000 1.0037 1.0073 1.0107 1.0136 1.0160 1.0174 1.0175

1.0000 1.0023 1.0045 1.0066 1.0084 1.0100 1.0110 1.0113

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

Columns 9 through 11 1.0000 1.0000 1.0000 1.0037 1.0023 1.0000 1.0073 1.0045 1.0000 1.0107 1.0066 1.0000 1.0136 1.0084 1.0000 1.0160 1.0100 1.0000 1.0174 1.0110 1.0000 1.0175 1.0113 1.0000 1.0155 1.0103 1.0000

1.0103 1.0072 1.0000 1.0000 1.0000 1.0000

>> contourf (uu, 'DisplayName', 'uu'); figure(gcf)

11109876543211234567891011

图一 超松弛

二 用块Jacobi迭代法求解线性方程组Au=f。 1. 基本原理：

对A做自然分解A=D-L-U=D-(L+U)

其中D是有A的对角线元素组成的矩阵，L是由A的对角线以下元素组成的矩阵，U是由A得对角线以上元素组成的矩阵。

于是将M=D,N=L+U,代入得到 Dx=(L+U)x+b

任取x的初值进行迭代

2. 算法：

（1）Gauss-Sediel迭代法原理

i?1 (k?1)(k?1)

xi?(bi?n?aj?1xjij???uaijxj)/aii,(k)j?i?1i?1,2,?,n五点差分格式：

(k?1)(k?1)(k?1)(k)(k)i,ji?1,ji,j?1i,j?1i?1,j i,j因为A可以写成块状，即：

u?(b?u?u?u)/4?A212??IA???????IA33????I????I??AN?1,N?1??T 如果把每一条线上的节点看作一个组v3?(u2,3,u3,3,...,uN,3)，可以把Au=f表示成块状求解：

Ajjvj?vj?1?vj?1?bj,

（2）计算步骤：

2?j?N 第一步：给场值u和场源b赋初值，及定义 第二步：用公式

Ajjvj?vj?1?vj?1?bjAii

，进行迭代计算

第三步：把第k次的u赋给ub，即ub=u；然后把第k+1次的u和ub进行比较，看是否达到精度，如果达到精度，则输出迭代次数k和精确解u。

3. 程序

[k,u]=kuai_GaussSeidel(n)

%用块Gauss-Sediel迭代法求解正方形域上的Poisson方程边值问题 %输入n,对x、y轴进行n等分；先确定场u的边界及场源b； %用k和u分别存放迭代次数和精确解

function [k,u]=kuai_GaussSeidel(n) %对x、y轴进行n等分 h=1/n; %步长

u=zeros(n+1,n+1); %对u赋初值 u(1,1:n+1)=1;u(n+1,1:n+1)=1;u(1:n+1,1)=1;u(1:n+1,n+1)=1;

b=zeros(n+1,n+1); %计算场源b for i=2:n+1 for j=2:n+1

b(i,j)=(i-1)*(j-1)*h^2; end end

b=h^2*b; for i=2:n

b(2,i)=b(2,i)+u(1,i);b(i,n)=b(i,n)+u(i,n+1);

b(n,i)= b(n,i)+u(n+1,i);b(i,2)= b(i,2)+u(i,1);

end

A=zeros(n-1,n-1); %定义矩阵的子块A for i=1:n-1

if i>1, A(i,i-1)=-1; end if i

e=0.000000001; %e是精度 tic;

for k=1:1000 %k是迭代次数 error=0; %误差

for j=2:n %进行迭代循环

ub=u; %ub是第k-1次迭代结果，用于和第k次迭代结果比较 if j==2, u(2:n,2)=pinv(A)*(u(2:n,3)+b(2:n,2)); end if j==n, u(2:n,n)=pinv(A)*(u(2:n,n-1)+b(2:n,n)); end

if j>2&j

error=max(max(abs(u-ub))); %error是前后两次迭代结果的对应元素的最大误差 if error<=e, break; end %判断误差是否达到精度 end

t=toc %统计程序运算时间

4. 计算结果： >> format short >> n=10;

>> [k,u]=kuai_GaussSeidel(n) t = 1.3280 k =

93 u =

Columns 1 through 8

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

1.0000 1.0011 1.0022 1.0031 1.0039 1.0044 1.0047 1.0045

1.0000 1.0022 1.0042 1.0061 1.0076 1.0087 1.0091 1.0088

1.0000 1.0031 1.0061 1.0088 1.0110 1.0126 1.0133 1.0128

1.0000 1.0039 1.0076 1.0110 1.0138 1.0159 1.0168 1.0162

1.0000 1.0044 1.0087 1.0126 1.0159 1.0189

1.0000 1.0047 1.0091 1.0133 1.0168 1.0203

1.0000 1.0045 1.0088 1.0128 1.0162 1.0201

1.0000 1.0037 1.0073 1.0107 1.0136 1.0175

1.0000 1.0023 1.0045 1.0066 1.0084 1.0113

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

Columns 9 through 11

1.0000 1.0000 1.0000 1.0037 1.0023 1.0000 1.0073 1.0045 1.0000 1.0107 1.0066 1.0000 1.0136 1.0084 1.0000 1.0160 1.0100 1.0000 1.0174 1.0110 1.0000 1.0175 1.0113 1.0000 1.0155 1.0103 1.0000 1.0103 1.0072 1.0000 1.0000 1.0000 1.0000

>> contourf (u, 'DisplayName', 'u'); figure(gcf)

1.0183 1.0194 1.0194 1.0208 1.0189 1.0203 1.0160 1.0174 1.0100 1.0110 1.0000 1.0000

11109876543211234567891011

图二 块Gauss-Sediel迭代法

三 （预条件）共轭斜量法求解线性方程组Au=f。 1.基本原理

（1）预条件共轭斜量法原理 ?1?TT?1Au?b?CACCu?Cb?Au?b

?1?T?1T其中A?CAC,b?Cb,u?Cu，

令M=CC称为预优矩阵，当M接近A时，A接近单位阵，cond(A)2接近1，对Au?f用共轭斜量法求解，可达到加速的目的。T（2）预优矩阵的选取

?A212??IA???????IA33????I????I??AN?1,N?1??

M取为块Jacobi迭代的块对角阵，?A22?M???????????An?1,n?1??A33?

2. 算法

(1)取初值u(0)?R(n),r(0)(0)?b?Au(0)(0),(0)(2)解方程组Mz(3)k?0,1,...,(0)?r,求出z，令p?z(0)?k?ur(k?1)(k?1)(z(k),r(k)))(k)(k)(k?1)(Ap?u?r(k)(k),p(k)??kp(k)??kAp(k?1)(k?1)(k)解方程组Mz?r)?k?(4)直到r(k?1)(z(k?1)(k),r,r(z)(k?1)p。(k?1)?z(k?1)??kp(k)??,输出u3.程序

%用J-PCG求解正方形域上的Poisson方程边值问题

%输入n,对x、y轴进行n等分；先确定场u的边界及场源b； %用k和u分别返回迭代次数和精确解 function [k,u]=J_CG(n)

h=1/n; %h步长

u=zeros(n+1,n+1); %对u赋初值

u(1,1:n+1)=1;u(n+1,1:n+1)=1;u(1:n+1,1)=1;u(1:n+1,n+1)=1; b=zeros(n+1,n+1); %计算场源b for i=2:n+1 for j=2:n+1

b(i,j)=(i-1)*(j-1)*h^2; end end b=h^2*b; for i=2:n

b(2,i)=b(2,i)+u(1,i);b(i,n)=b(i,n)+u(i,n+1);

b(n,i)= b(n,i)+u(n+1,i);b(i,2)= b(i,2)+u(i,1);

end

z=zeros(n-1,n-1);r=zeros(n-1,n+1);p=zeros(n-1,n-1);bb=zeros(n-1,n-1); A=zeros(n-1,n-1); %定义矩阵的子块A for i=1:n-1

if i>1, A(i,i-1)=-1; end if i

end

for j=2:n

if j==2, bb(:,1)=A*u(2:n,2)-u(2:n,3); end ?=A*u

if j==n, bb(:,n-1)=A*u(2:n,n)-u(2:n,n-1); end

if j>2&j

r=b(2:n,2:n)-bb; %计算r0 for i=1:n-1 %计算z0 z(:,i)=pinv(A)*r(:,i); end

p=z; %计算p0 tic;

e=0.000000001; %e是精度 for k=1:1000

error=0;a=0;aa=0;cc=0;ub=u; for i=1:n-1

aa=aa+dot(z(:,i),r(:,i));

if i==1, cc=cc+dot((A*p(:,i)-p(:,i+1)),p(:,i));end

if i>1&i

a=aa/cc; %a是(z(k),r(k))/(Ap(k),p(k))，确定搜索步长 u(2:n,2:n)= u(2:n,2:n)+a*p; %对u进行更新计算

rr=r;zz=z; %rr和zz存储第k-1次的r和z for i=1:n-1

if i==1, r(:,i)= r(:,i)-a*(A*p(:,i)-p(:,i+1));end

if i>1&i

z(:,1:n-1)=pinv(A)*r(:,1:n-1); %Mz kk=0;mm=0;

for i=1:n-1

kk=kk+dot(z(:,i),r(:,i)); mm=mm+dot(zz(:,i),rr(:,i)); end

beita=kk/mm; ?ita是(z(k?1)(k?1)?r(k?1)

,r(k?1))/(z(k),r(k))

p=z+beita*p; %对p进行更新，确定搜索方向 error=sqrt(norm(u-ub)); %error是统计误差

if error<=e, break; end %判断误差是否达到精度 end

t=toc %统计计算时间

4. 计算结果： >> format short >> n=10;

>> [k,u]=J_CG(n) t = 0.1100 k = 37 u =

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

1.0000 1.0011 1.0022 1.0031 1.0039 1.0045 1.0037 1.0023 1.0000

1.0000 1.0022 1.0042 1.0061 1.0076 1.0088 1.0073 1.0045 1.0000

1.0000 1.0031 1.0061 1.0088 1.0110 1.0128 1.0107 1.0066 1.0000

1.0000 1.0039 1.0076 1.0110 1.0138 1.0162 1.0136 1.0084 1.0000

1.0000 1.0044 1.0087 1.0126 1.0159 1.0189 1.0160 1.0100 1.0000

1.0000 1.0047 1.0091 1.0133 1.0168 1.0203 1.0174 1.0110 1.0000

1.0000 1.0045 1.0088 1.0128 1.0162 1.0201 1.0175 1.0113 1.0000

1.0000 1.0037 1.0073 1.0107 1.0136 1.0175 1.0155 1.0103 1.0000

1.0000 1.0023 1.0045 1.0066 1.0084 1.0113 1.0103 1.0072 1.0000

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

>> contourf (u, 'DisplayName', 'u'); figure(gcf)

1.0000 1.0000 1.0044 1.0047 1.0087 1.0091 1.0126 1.0133 1.0159 1.0168 1.0183 1.0194 1.0194 1.0208 1.0189 1.0203 1.0160 1.0174 1.0100 1.0110 1.0000 1.0000

11109876543211234567891011

图三 J-PCG

四 三种迭代法效率分析

有场的等值图可以看出，三种迭代方法的结果（达到相同的精度）比较一致，但是J-PCG只需要37次（耗时0.1100秒）迭代即可达到比较好的结果；超松弛迭代法需要48次（耗时0.0310秒）迭代即可达到满意的结果；块Gauss-Sediel迭代法需要93次（耗时1.3280）迭代才到达满意的结果。所以对此问题来说，超松弛迭代法和J-PCG法的效率要好于块Gauss-Sediel迭代法。

一. 任务：

用MATLAB语言编写连续函数最佳平方逼近的算法程序（函数式M文件）。并用此程

序进行数值试验，写出计算实习报告。 二. 程序功能要求：

在书中Page355或Page345的程序leastp.m（见附一）的基础上进行修改，使其更加完善。要求算法程序可以适应不同的具体函数，具有一定的通用性。所编程序具有以下功能：

1. 用Lengendre多项式做基，并适合于构造任意次数的最佳平方逼近多项式。

P0(x)?1,P1(x)?x可利用递推关系 Pn(x)???(2n?1)xPn?1(x)?(n?1)Pn?2(x)??/n

n?2,3,.....2. 被逼近函数f(x)不用内联函数构造，而改用M文件建立数学函数。这样，此程序可通过

修改建立数学函数的M文件以适用不同的被逼近函数（要学会用函数句柄）。 3. 要考虑一般的情况f(x)?[a,b]?[?1,?1]。因此，程序中要有变量代换的功能。 4. 计算组合系数时，计算函数的积分采用变步长复化梯形求积法（见附三）。

5. 程序中应包括帮助文本和必要的注释语句。另外，程序中也要有必要的反馈信息。 6. 程序输入：（1）待求的被逼近函数值的数据点x0(可以是一个数值或向量)

（2）区间端点：a,b。

7. 程序输出：（1）拟合系数：c0,c1,c2,...,cn

（2）待求的被逼近函数值

1. 基本原理：

设内积空间Ｘ，取线性无关组?span{?0?x?，?，?1?x?，?n?x?,构造Ｘ的子空间 Ｈ＝

0?x?，?，?1?x?，?n?x?,},对被逼近函数f(x)不属于Ｈ构造逼近函数s(x):

ns(x)= ?j?02

c??x?

jj求出系数c0,c1, …,cn,使‖f(x)-s(x)‖2=min,则函数s(x)就是f(x)H在中的最佳平方逼近元。

若取正交函数组?0?x?，?，?1?x?，?n?x?,对权函数ρ(x)满足

?x?,?j?x?=???x???x???x?dx=?ijb（i,j=0,1,…，n）

??i?b0i?j2?i?x?2i?j

以Legendre正交多项式为函数的最佳平方逼近多项式，若f(x)=xcos(x),x∈[-1，1]，用p1(x)=1,p2(x)=x,p3(x)=(3x2-1)/2,p4(x)=(5x3-3x)/2构造函数f(x)的三次最佳平方逼近多项式 :

s(x)=c1p1(x)+c2 p2(x)+c3 p3(x)+c4 p4(x)

2 程序功能：

该程序可对被拟合函数做以Legendre多项式为基的最佳平方拟合，拟合数据可以是一个点或一个向量，并绘出该函数与拟合多项式以及泰勒展式在指定区间上的对比图。 3 使用说明：

程序用法简单，只需输入拟合数据点x，拟合区间a、b的值，以及[c,s]=leasp(x,a,b),程序即可自动运行。 4 源程序： (1) 主程序M文件

%LEASTP.M: 以Legendre多项式为基的最佳平方拟合 function [c,s]=leastp(x,a,b) %x:输入拟合点 %func:输入拟合函数 %c:输出拟合系数 %s:输出拟合值

if nargin<3, disp('必须输入三个参数'); return;end ff=@func; for i=1:4

pp(i)=2/(2*(i-1)+1); end

c(1)=quad('func',a,b)/pp(1)*(2/(b-a));

c(2)=quad('fp2',a,b,[ ],[ ],a,b)/pp(2)*(2/(b-a)); c(3)=quad('fp3',a,b,[ ],[ ],a,b)/pp(3)*(2/(b-a)); c(4)=quad('fp4',a,b,[ ],[ ],a,b)/pp(4)*(2/(b-a)); s=c(1)+c(2)*p2(x,a,b)+c(3)*p3(x,a,b)+c(4)*p4(x,a,b); x0=a:0.15:b;

s0=c(1)+c(2)*p2(x0,a,b)+c(3)*p3(x0,a,b)+c(4)*p4(x0,a,b);

plot(x0,s0,'r-'); %绘制多项式拟合图 hold on

fplot(ff,[a,b,feval(ff,a),feval(ff,b)],'c'); % 绘制拟合函数图 syms xx

ff1=xx*cos(xx); y1=taylor(ff1,xx,4); hold on

ezplot(y1,[a,b]); % 绘制拟合函数的泰勒展式图像 title('comparison figure of fitted funcyion,fitting polynomials and taylor formulation');

plot(x,s,'r *'); %数据点

legend('原 函 数 f(x)','最佳逼近s(x)','泰勒展开T(x)','数据点');

(2) legendre多项式M文件 function y=p2(x,a,b) t=(2*x-a-b)./(b-a); y=t;

function y=p3(x,a,b) t=(2*x-a-b)/(b-a); y=(3*t.^2-1)/2; function y=p4(x,a,b) t=(2*x-a-b)./(b-a); y=(5*t.^3-3*t)/2;

(3) 函数与多项式内积M文件 function y=fp2(x,a,b)

y=feval('p2',x,a,b).*feval('func',x); function y=fp3(x,a,b)

y=feval('p3',x,a,b).*feval('func',x); function y=fp4(x,a,b)

y=feval('p4',x,a,b).*feval('func',x);

(4) 函数M文件 function f=func(x) f=x.*cos(x);

本文来源：https://www.bwwdw.com/article/odg3.html