大学生科技创新论文题目：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:nA(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 00.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:nd(k,j)=(d(k,j-1)-d(k-1,j-1))/(x(k)-x(k-j+1)); end endc=d(n,n);for k=(n-1):-1:1c=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）对上述数据做线性拟合，关键在于找出法方程。

定义如下的内积函数：060(,)i i iw ϕϕ==∑，061(,)i i i iw x ϕϕ==∑，6211(,)i ii iw x ϕϕ==∑，61(,)f i i i iw 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 enda f fb 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 00.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 054320.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 =和11D 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 = 28375H =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=-LB=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 0B = 80707/5C = 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 1x1=B*x0+f k=k+1 fprintf('%2d ',k);disp(x1);if norm(x1-x0,inf)<eps;fprintf('\n'); 第次迭代的近似解满足精度要求的近似解 disp(x1); break end if k>Nfprintf('\n'); break endx0=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 1x1=B*x0+f k=k+1 fprintf('%2d ',k);disp(x1);if norm(x1-x0,inf)<eps; fprintf('\n');disp(x1);break end 第次迭代的近似解满足精度要求的近似解if k>Nfprintf('\n');break endx0=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 endU=-(H-D-L);L=-L;超松弛迭代2while 1k=k+1;x1=B*x0+f;if x1-x0<1e-5break;else((()*()));2/(1(1));(*)*((1)**);*(*)*; x0=x1;endendF=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.3473k = 3219x1 =-3.999859.9974-179.9936139.9958可以看出，条件数小于上述迭代法，结论成立。