大学生科技创新论文

大学生科技创新论文

题目：MATLAB与计算方法关于数值计算问题姓名：张欣（1008300062）杜昕阳（1008300061）

班级：10级信息班

专业：信息与计算科学

内蒙古包头师范学院

二零一二年十一月一日

一、题目：

MATLAB与计算方法关于数值计算问题

二、摘要：

MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。[1] MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连matlab开发工作界面

接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++，JA V A的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。

《计算方法》比较全面地介绍了科学与工程计算中常用的计算方法，具体介绍了这些计算方法的基本理论与实际应用，同时对这些数值计算方法的计算效果、稳定性、收敛效果、适用范围以及优劣性与特点也作了简要的分析。全书共9章，内容包括引论、线性代数方程组求解方法、非线性方程求根、函数插值、函数逼近、矩阵特征值与特征向量的数值算法、数值积分与数值微分、常微分方程初值问题的数值解法、自治微分方程稳定区域的计算等。《计算方法》概念清晰，语言叙述通俗易懂，理论分析严谨，结构编排由浅入深，在分析问题时注重启发性，例题选择具有针对性且注重实际应用。前8章附有一定数量的习题，供读者学习时进行练习。《计算方法》可作为高等院校数学与应用数学、信息与计算科学、应用物理学、计算机科学等专业的高年级本科生和工科硕士研究生使用，也可供从事科学与工程计算的科技工作者参考。

三、引言：

计算方法又称“数值分析”。是为各种数学问题的数值解答研究提供最有效的算法。主要内容为函数逼近论，数值微分，数值积分，误差分析等。常用方法有迭代法、差分法、插值法、有限元素法等。现代的计算方法还要求适应电子计算机的特点。数值分析即“计算方法”.

四、 正文：

插值与拟合 已知()101x f x e x =+-

a) 求函数在0,0.2,0.4,0.6,0.8,1x =处的函数值；（先编函数，再求值保存到

向量中）

b) 对上述数据进行多项式插值，作出多项式5()y P x =的图像，与原函数图

象比较；

（先列出差分表，再用牛顿插值公式编写出多项式函数） c) 对上述数据做线性拟合，作出多项式1()y P x =的图像；

（先定义内积函数，再列出法方程，然后求解，最后编出多项式函数）

d) 构造()f x 在[0,1]区间内的5次切比雪夫多项式5()T x ，并作出图像。

（先生成切比雪夫点，再列出差分表，再插值）

（1）函数代码为 %function y=P(x)y=exp(x)+10*x-1;函数文件 %x=[0:0.2:1];y=P(x)相应点的函数值

结果为 y = 0 2.2214 4.4918 6.8221 9.2255 11.7183

（2）差分表代码为 %function A=ff(x,y)n=length(x);A=zeros(n,n+1);A(:,1)=x';

A(:,2)=y';for j=3:n for i=j-1:n

A(i,j)=A(i,j-1)-A(i-1,j-1); end end 差分表代码 x=[0:0.2:1];

y=P(x);A=ff(x,y)

结果为 A =

0 0 0 0 0 0 0 0.2000 2.2214 2.2214 0 0 0 0

0.4000 4.4918 2.2704 0.0490 0 0 0 0.6000 6.8221 2.3303 0.0599 0.0109 0 0 0.8000 9.2255 2.4034 0.0731 0.0133 0.0024 0 1.0000 11.7183 2.4927 0.0893 0.0162 0.0029 0

牛顿插值公式代码为 %x=0:0.2:1;y=P(x);C=newpoly(x,y);

xx=0:0.1:1;

yy=polyval(C,xx);yyy=poly2str(C,'x')plot(x,y,'ro',xx,yy)

牛顿插值及相应图像

%function c=newpoly(x,y)n=length(x);d=zeros(n,n);d(:,1)=y';

for j=2:n for k=j:n

d(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1)); end end

c=d(n,n);

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

c=conv(c,poly(x(k))); m 牛顿插值法程序

=length(c);

c(m)=c(m)+d(k,k);end

结果为 54320.0138540.0348660.170410.4990711.0001yyy x x x x x =++++

（3）对上述数据做线性拟合，关键在于找出法方程。定义如下的内积函数：

06

0(,)i i i

w ??==∑，06

1(,)i i i i

w x ??==∑，6

211(,)i i

i i

w x ??==∑，6

1(,)f i i i i

w f ??==∑，这里取

权重为1，然后再构造如下的法方程，1

(,)k j j k j a d ??==∑（k=0,1）,由法方程得到

线性方程组，从而通过解线性方程组得到线性拟合曲线的系数。

线性拟合程序代码为 (1)(,,)();(1,1);

(,1)c=a\b;

c=fli ;

1:1

(:pud(,).;*;

*;c);k function c Ispoly x y m n length x b zeros m f zeros n m for k m f k x end

a f f

b f y -===+=+=+'='=''= x=0:0.2:1;y=P(x);

C=Ispoly(x,y,3);

xx=0:0.1:1;

yy=polyval(C,xx);

yyy=poly2str(C,'x')

plot(x,y,'ro',xx,yy)

结果为 320.279240.424311.01470.00018852yyy x x x =++-

（4）利用6()T x 的零点作插值节点，构造()f x 在[0,1]区间内的5次多项式。由于切比雪夫零点是在区间[1,1]-之间，对于一般区间只要利用如下的变换即可。

()/2*cos((2*1)*/2*(1))()/2,0,1....k x b a k pi n a b k n =-++++= 所

以首先要生成切比雪夫点。，利用上述区间变换的公式可得节点：

0124560.98,0.85,0.63,0.37,0.15,0.017x x x x x x ======，进一步得到差分表，找

出插值函数，并作出相关图像。

代码为 x=[0.98,0.85,0.63,0.37,0.15,0.017];y=P(x);A=ff(x,y)C=newpoly(x,y);xx=0:0.1:1;yy=polyval(C,xx);

yyy=poly2str(C,'x')plot(x,y,'ro',xx,yy)

，

结果为A =

0.9800 11.4645 0 0 0 0 0 0.8500 9.8396 -1.6248 0 0 0 0

0.6300 7.1776 -2.6620 -1.0372 0 0 0 0.3700 4.1477 -3.0299 -0.3678 0.6694 0 0 0.1500 1.6618 -2.4859 0.5440 0.9118 0.2424 0 0.0170 0.1871 -1.4747 1.0112 0.4672 -0.4446 0

54320.0138540.0349150.170310.4991211.0001 1.0643006yyy x x x x x e =++++--

4、线性方程组的迭代解法

（1）构造4阶Hilbert 矩阵H ，求()cond H ∞；

（2）分裂H D L U =--，求1()cond D H -∞，1(())cond D L H -∞-；

（3）设T [1,1,1]b =，找ω使得迭代法(1)()()()k k k x x b Hx ω+=+-收敛速度最快； （4）分别用（3）中的方法，雅可比迭代法，高斯塞德尔迭代法，使

(1)()15k k x x e +∞

-<- 时终止，比较迭代次数与（1）

（2）题中条件数的关系。

（5）用共轭梯度法分别求解Hx b =和11

D Hx D b --=，(1)()

15k k x x e +∞

-<-终

止，比较迭代次数与系数矩阵条件数的关系。

（1）代码为 %Hilbert format rat H=hilb(4)四阶矩阵 %Hilbert H=hilb(4);A=norm(H,inf)*norm(inv(H),inf)四阶矩阵的无穷范数

结果为 A = 28375

H =

1 1/

2 1/

3 1/

4 1/2 1/3 1/4 1/

5 1/3 1/4 1/5 1/

6 1/4 1/5 1/6 1/7

（2）代码为

%H D-L-U H=hilb(4);D=diag(diag(H))L=tril(H);for i=1:4 for j=1:4 if i==j L(i,j)=0; end end end U=-(H-D-L)L=-L

B=norm(inv(D)*H,inf)*norm(inv(inv(D)*H),inf)C=norm(inv(D 分裂为，并求范数-L)*H,inf)*norm(inv(inv(D-L)*H),inf) 结果为

D =

1 0 0 0 0 1/3 0 0 0 0 1/5 0 0 0 0 1/7 U =

0 -1/2 -1/3 -1/4 0 0 -1/4 -1/5 0 0 0 -1/6 0 0 0 0 L =

0 0 0 0

-1/2 0 0 0 -1/3 -1/4 0 0 -1/4 -1/5 -1/6 0

B = 80707/5

C = 8225/3

（3）要使迭代速度最快，只要有w 满足迭代矩阵的谱半径最小，则代码为：

function y=f(w)H=hilb(4);

y=max(abs(eig(eye(4)-w*H))); w=fminsearch('f',0)

结果为w = 1.3330

（4）雅阁比迭代法和高斯赛德尔迭代法程序为

%A=hilb(4)';b=[1,1,1,1]';

if(any(diag(A))==0);error('')end eps=input('eps=');N=input('N=')D=diag(diag(A)); B=inv(D)*(D-A); f=inv(D)*b;k=0; x0=zeros(size(b)); 雅阁比迭代主对角线上存在零元，迭代程序终止请输入误差限请输入最大迭代次数while 1

x1=B*x0+f k=k+1 fprintf('%2d ',k);disp(x1);

if norm(x1-x0,inf)

fprintf('\n'); 第次迭代的近似解满足精度要求的近似解 disp(x1); break end if k>N

fprintf('\n'); break end

x0=x1; end

迭代次数超限

%A=hilb(4)'; b=[1,1,1,1]';

if(any(diag(A))==0);error('')end eps=input('eps=');N=input('N=')

D_L= tril(A);U=-(A-D_L); B=inv(D_L)*U;f=inv(D_L)*b;k=0;x0=zeros(size(b));高斯赛德尔迭代主对角线上存在零元，迭代程序终止请输入误差限请输入最大迭代次数 while 1

x1=B*x0+f k=k+1 fprintf('%2d ',k);disp(x1);

if norm(x1-x0,inf)N

fprintf('\n');break end

x0=x1; end

迭代次数超限

代码为

%function [k,x1,F]=SOR(x0)H=hilb(4);b=[1,1,1,1]';D=diag(diag(H));L=tril(H);for i=1:4 for j=1:4 if i==j L(i,j)=0; end end end

U=-(H-D-L);L=-L;

超松弛迭代

2while 1 k=k+1; x1=B*x0+f; if x1-x0<1e-5 break; else

((()*()));2/(1(1));

(*)*((1)**);*(*)*; x0=x1; end end

F=norm(B,in 0;f)*p max eig inv D L U w sqrt p B inv D w L w D w U f w inv D w L b k =+=+-=--+=-=norm(inv(B),inf)

[k,x1]=SOR([1,1,1,1]')

结果为 F = 76.3473

k = 3219x1 = -3.9998 59.9974 -179.9936

139.9958

可以看出，条件数小于上述迭代法，结论成立。

（5）共轭梯度法是构造一组正交向量，使得搜索方向快速下降，对于最速下降法的一种改进方法。 程序代码为

%function [k,x1]=CG(A,b,x0)r=b-A*x0;p=r;AP=A*p;k=0;while k<1000

k=k+1;a=r'*r/(p'*AP);x1=x0+a*p;r1=r-a*AP;B=r1'*r1/(r'*r);p1=r1+B*p;F=norm(x1-x0,inf);if F<1e-5break;end r=r1;x0=x1;p=p1;AP=A*p;en 共轭梯度法算法程序d

b=[1,1,1,1]';

x0=[0,0,0,0]';

A=hilb(4);

[k,x1]=CG(A,b,x0)

b=[1,1,1,1]';

x0=[0,0,0,0]';

H=hilb(4);

D=diag(diag(H));

A=inv(D)*H;

A=A'*A;

b=inv(D)*b;

b=b'*b;

[k,x1]=CG(A,b,x0)

T1=norm(A,inf)*norm(inv(A),inf) T2=norm(H,inf)*norm(inv(H),inf)

结果为

k = 5

x1 =

-4.0000 60.0000 -180.0000 140.0000k = 6

x1 = 1.0e+006 * -0.2513

2.7485

-6.5197

4.2000

条件数为T1 =1.5349e+008 T2 = 2.8375e+004

五、结论：

由于该题所给数据局限，两种迭代方法无法迭代出满意的结果，尤其是雅阁比迭代法，竟然是不收敛的，但是可以知道一个结论，就是条件数越小，迭代速度越快。所以，如果需要开发或者衡量一种迭代方法的话，可以考虑相应的条件数，令其足够小。而第三种迭代法，条件数小于上述迭代法，且迭代结果较为合理。

六、参考文献：

1.李庆阳、王能超、易大义.数值分析，第三版。武汉：华中理工大学出版社，

1986

2.冯康等编，数值计算方法，北京：国防工业出版社，1978

3.黄友谦，李岳生，数值逼近，第二版，北京：高等教育出版社，1987

4.J.Stoer,R.Bulirsch.Introduction to Numerical Analysis ,Spring-Verlay, New

York,1980

5.De Boor carl.A Practical Guide to Spline .Spring-Verlag,New York ,1978

6.J.M奥特加，数值分析，张丽君等译，北京：高等教育出版社，1973

7.菜大用，数值代数，北京：清华大学出版社，1987

8.路易斯.汉格曼,大卫.杨，实用迭代法.菜大用等译，北京：清华大学出版社，

1981

9.易大义、李江云.插值、逼近、曲线和曲面拟合，上海：复旦大学出版社，1992

10.A.R，高尔腊伊，G.A，瓦特桑.矩阵特征值问题的计算方法，唐换文等译，上

海：上海

科技出版社，1980

11.李庆阳、莫子忠、祁立群，非线性方程组的数值解法，北京：科学出版社，

1987