《数值分析matlab程序.doc》由会员分享，可在线阅读，更多相关《数值分析matlab程序.doc（46页珍藏版）》请在三一办公上搜索。
1、数值分析（matlab程序）曹德欣 曹璎珞第一章 绪论数值稳定性程序，计算P11 试验题一积分function try_stable global n global a a=input(a=); N = 20; I0 = log(1+a)-log(a); I = zeros(N,1); I(1) = -a*I0+1; for k = 2:N I(k) = - a*I(k-1)+1/k; end II = zeros(N,1); if a=N/(N+1) II(N) = (1+2*a)/(2*a*(a+1)*(N+1); else II(N) =(1/(a+1)*(N+1)+1/N)/2; en
2、d for k = N:-1:2 II(k-1) = ( - II(k) +1/k) / a; end III = zeros(N,1); for k = 1:N n = k; III(k) = quadl(f,0,1); end clc fprintf(n 算法1结果 算法2结果 精确值) for k = 1:N, fprintf(nI(%2.0f) %17.7f %17.7f %17.7f,k,I(k),II(k),III(k) endfunction y = f(x) global n global a y = x.n./(a+x); return第二章 非线性方程求解下面均以方程y=x
3、4+2*x2-x-3为例：1、二分法function y=erfen(a,b,esp) format long if nargin3 esp=1.0e-4; end if fun(a)*fun(b)esp if fun(a)*fun(c)0 b=c; c=(a+b)/2; elseif fun(c)*fun(b)=1.0e-4) & (n=1.0e-4) & (n=100000000) %判断两个条件截止 x0=x1; %将x1赋给x0 x1=x2; %将x2赋给x1 x2=x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0); %迭代运算 n=n+1; end y=x2 n
4、function y=fun(x) y=x4+2*x2-x-3;第四章题目：推导外推样条公式：，并编写程序与Matlab的Spline函数结果进行对比，最后调用追赶法解方程组。解：设的二阶导数值，因为在上是三次多项式，所以在上是一次多项式，且可表示为 （1）对式（1）进行两次积分可得： （2）由式（2）根据和可得系数为：， （3）将式（3）代入式（2）可得为：分别对和进行求导，并根据两者相等整理得： （4）其中，但是要解给定的方程组，还需要两个另外的条件，而外推样条插值的条件可通过下面推理得出：令式（1）中的等于，可得： （5） （6）令式（4）中，然后将式（5）代入得： （7）令式（4）中，
5、然后将式（6）代入得： （8）将方程组（4）和式（7）、式（8）联立展开即得题目所求。按照推导出的外推样条插值公式编程可得如下M文件（spline_wt.m）function spline_wtX=0 1 2 3 4 5;Y=0 0.5 2 1.5 3.5 1.9;% 调用自编程序S = spline_w(X,Y); %调用matlab提供程序S1=spline(X,Y);fnplt(S,r,2); % 作图hold onfnplt(S1,b,1);hold onplot(X,Y,or); % 画上节点title(红线为自编程序曲线，蓝线为自带程序曲线)% pp的第j行表示第j个三次多项式的4
6、系数并写出分段多项式pp = S.coefs P1 = poly2str(pp(1,:),(x-0) P2 = poly2str(pp(2,:),(x-1) P3 = poly2str(pp(3,:),(x-2)ppp=S1.coefs PP1 = poly2str(ppp(1,:),(x-0) PP2 = poly2str(ppp(2,:),(x-1) PP3 = poly2str(ppp(3,:),(x-2)function sp = spline_w(X,Y)n = length(X); h = diff(X); d = diff(Y)./h; d1(2:n-1)=6*diff(d)./
7、(h(1:n-2)+h(2:n-1); mu(2:n-1)=h(1:n-2)./(h(1:n-2)+h(2:n-1); mu(n-1)=1-h(n-2)/h(n-1); la(2:n-1)=1-mu(2:n-1); la(2)=1-h(1)/h(2);% 计算三对角方程组 a = mu(3:n-1); b = 2*ones(n-2,1); b(1)=2+h(1)/h(2); b(n-2)=2+h(n-2)/h(n-1); c = la(2:n-2); u(1:n-2) = d1(2:n-1);% 调用追赶法解方程组 M(2:n-1) = tridiag(a,b,c,u); M(1)=(1+h(
8、1)/h(2)*M(2)-M(3)*h(1)/h(2);M(n)=(1+h(n-1)/h(n-2)*M(n-1)-M(n-2)*h(n-1)/h(n-2);% 下面计算分段多项式的四个系数S=zeros(n-1,4);for k=0:n-2S(k+1,1)=(M(k+2)-M(k+1)/(6*h(k+1);S(k+1,2)=M(k+1)/2;S(k+1,3)=d(k+1)-h(k+1)*(2*M(k+1)+M(k+2)/6;S(k+1,4)=Y(k+1);end sp = mkpp(X,S);%转换成 Matlab 格式% - 求解三对角线性方程组的追赶法 -function x=tridia
9、g(a,b,c,d)n = length(d); x = zeros(n,1); for k = 2:n b(k) = b(k) - a(k-1)*c(k-1)/b(k-1); d(k) = d(k) - a(k-1)*d(k-1)/b(k-1); end x(n) = d(n)/b(n); for k = n-1:-1:1 x(k) = ( d(k)-c(k)*x(k+1) ) / b(k);end根据所编写的程序可得三次多项式及系数为（包括自编和提供）：pp = -0.9511 3.3533 -1.9022 0 -0.9511 0.5000 1.9511 0.5000 1.7556 -2.
10、3533 0.0978 2.0000 -1.5711 2.9133 0.6578 1.5000 -1.5711 -1.8000 1.7711 3.5000P1 =-0.95111 (x-0)3 + 3.3533 (x-0)2 - 1.9022 (x-0)P2 =-0.95111 (x-1)3 + 0.5 (x-1)2 + 1.9511 (x-1) + 0.5P3 =1.7556 (x-2)3 - 2.3533 (x-2)2 + 0.097778 (x-2) + 2ppp = -0.9511 3.3533 -1.9022 0 -0.9511 0.5000 1.9511 0.5000 1.7556
11、 -2.3533 0.0978 2.0000 -1.5711 2.9133 0.6578 1.5000 -1.5711 -1.8000 1.7711 3.5000PP1 =-0.95111 (x-0)3 + 3.3533 (x-0)2 - 1.9022 (x-0)PP2 =-0.95111 (x-1)3 + 0.5 (x-1)2 + 1.9511 (x-1) + 0.5PP3 =1.7556 (x-2)3 - 2.3533 (x-2)2 + 0.097778 (x-2) + 2根据所编写的程序可得生成的图形如图1所示第五章摘要：针对多项式检验的困难性，文中先将多项式线性化为多元线性问题。在此基
12、础上，简要推导出多元线性拟合的数学模型及系数的最小二乘估计（以二元为例），同时为了检验建立模型的准确性，推导出检验模型的三个指标（误差差方和、相关检验和F检验）计算公式，最后用某矿的煤样实例进行试验研究，得出各次多项式拟合的检验结果。1 多元线性拟合的数学模型设一个随机变量与个一般变量的内在联系是线性的，而且统计相关。与这些之间可用如下线性关系表示： （1）上式称为多元线性拟合的数学模型。式中是个待估参数，这些参数被称为回归系数；是服从的随机变量。多项式拟合可以通过线性化来化为多元线性拟合，如多项式模型为： 对上式线性化，令，即可得如式（1）所示的多元线性模型，因此我们只需讨论多元线性模型即可
13、。为了方便我们采用二次多项式来推导，多次多项式可根据二次近似得到。2 系数的最小二乘估计设给定组实测数据，将其带入式（1）可得个观测方程： （2）式中随机误差是个相互独立且服从的随机变量。令，则式（2）可写为：设是采用最小二乘估计方法求得的估值，则多元线性方程为：，则与实测值之间的差值为：用最小二乘估计，所有实验点上偏差平方和最小的要求下，设，则可得法方程用矩阵表示为：于是， （3）由（3）式解出，即可得多元线性方程模型。3方程模型的检验在实际问题中，事先并不能断定随机变量与一般变量之间是否有线性关系。用组实测数据按最小二乘的方法总能找出回归方程模型。但所得的方程模型是否具有实际意义，需要对所
14、得到的方程模型进行统计检验。本文采用了误差平方和和两种统计方法进行检验。3.1 误差平方和按下式可计算出总差方和：，自由度为因为：=0则可分解为和，称为回归差方和，反映自变量的重要程度自由度为；为误差差方和，反映试验误差及其他因素对试验结果的影响，自由度为。3.2 相关系数检验法一元线性方程中相关系数定义为：，其中，分别表示的均方差和协方差。因为它们是未知的，故用子样均方差和协方差来代替。最后得相关系数的估值： 其中，分别为子样的方差和协方差， 在多元线性方程中，定义复相关系数：，当时，线性关系密切，而当时，线性关系不密切，甚至不存在线性关系。可通过下式计算来得到值： 3.3 F检验如果变量与
15、之间无线性关系，则模型的一次项系数均为零。所以，要检验变量与变量之间是否有线性关系，就是要检验假设是否成立从，得定义可看出两者是独立的，且，则在矩阵满秩和原假设成立的条件下在水平下，若计算的（查F分布表），认为不可信，方程模型效果显著；反之，则认为可信，线性效果不显著，所求方程没有实际意义。4 MATLAB实现及分析实验中取用某煤矿的18个煤样，所得的18组数据列于表1，分析容重和灰分之间的关系。表1取样数据表样品号123456789101112容重x1.51.21.71.41.81.31.31.51.71.31.51.5灰分y%25430203675243341724样品号131415161
16、718容重x1.61.41.61.51.41.5灰分y%2562624209根据上述数据利用matlab的多项式进行拟合，并得到各次拟合的检验结果如表2，表2 多项式拟合检验结果表多项式次数1234临界值误差差方和378.2118373.5294369.4278358.566相关检验0.8930.89440.89560.89880.47F检验62.961129.97718.911213.66774.49拟合结果图形见图1图1 多项式拟合效果图从表2中可以看出，随着次数的增加，误差差方和逐渐减少，说明拟合越接近于实际曲线。这也可以通过相关检验的变化得以证明，次数越高，相关系数越接近于1且大于临界
17、值，说明多项式线性化后的模型显著且逐渐接近于线性关系。而F检验的结果（方程显著且大于临界值）却说明当次数增加到一定程度时，拟合模型可能出现不显著，也即拟合模型不可用。从上述分析可以看出，对于一个特定的问题，并不是建立模型的次数越多越接近于真实曲线。当多项式的次数超过F检验结果可用的情况下，高次数的多项式拟合反而成为一种错误的结果。因此，当我们在建立模型时，对建立的模型结果要进行检验，不能只根据某些实验数据拟合结果可用，就用它来进行下一步的预测研究。附录：matlab源码function poly_nihe%多元线性拟合及检验 x=1.5 1.2 1.7 1.4 1.8 1.3 1.3 1.5
18、1.7 1.3 1.5 1.5 1.6 1.4 1.6 1.5 1.4 1.5;y=25 4 30 20 36 7 5 24 33 4 17 24 25 6 26 24 20 9;nn=length(x);n=1; p=polyfit(x,y,n); xi=linspace(1.2,1.8,100); z=polyval(p,xi); subplot(2,2,1),plot(x,y,o,xi,z) title(一次多项式拟合) ay=mean(y); dy=y-ay; q=0; for i=1:nn q=q+dy(i)2; end ax=mean(x); dx=x-ax; qx=0; for
19、i=1:nn qx=qx+dx(i)2; end cxy=0; for i=1:nn cxy=cxy+dx(i)*dy(i); end r=cxy/sqrt(qx*q) %自由度为16及显著性水平为0.05查相关系数临界值表，得0.468xp=1.5 1.2 1.7 1.4 1.8 1.3 1.3 1.5 1.7 1.3 1.5 1.5 1.6 1.4 1.6 1.5 1.4 1.5;yp1=polyval(p,xp); v=y-yp1; q2=0; for i=1:nn q2=q2+v(i)2; end q2 q1=q-q2; f=q1*(nn-2)/q2 %自由度为16及显著性水平为0.0
20、5查F分布表，得4.49n=2; p=polyfit(x,y,n); xi=linspace(1.2,1.8,100); z=polyval(p,xi); subplot(2,2,2),plot(x,y,o,xi,z) title(二次多项式拟合) %F检验 yp2=polyval(p,xp); v=y-yp2; q2=0; for i=1:nn q2=q2+v(i)2; end q2 q1=q-q2; f=q1*(nn-n-1)/(q2*n) %R相关系数检验 r_2=q1/q; ff=r_2*(nn-n-1)/(1-r_2)*n); r=sqrt(n*f/(nn-n-1+n*f)n=3;
21、p=polyfit(x,y,n); xi=linspace(1.2,1.8,100); z=polyval(p,xi); subplot(2,2,3),plot(x,y,o,xi,z) title(三次多项式拟合) yp3=polyval(p,xp); v=y-yp3; q2=0; for i=1:nn q2=q2+v(i)2; end q2 q1=q-q2; f=q1*(nn-n-1)/(q2*n) %R相关系数检验 r_2=q1/q; ff=r_2*(nn-n-1)/(1-r_2)*n); r=sqrt(n*f/(nn-n-1+n*f)n=4; p=polyfit(x,y,n); xi=l
22、inspace(1.2,1.8,100); z=polyval(p,xi); subplot(2,2,4),plot(x,y,o,xi,z) title(四次多项式拟合) yp4=polyval(p,xp);v=y-yp4;q2=0; for i=1:nn q2=q2+v(i)2; end q2 q1=q-q2; f=q1*(nn-n-1)/(q2*n) %R相关系数检验 r_2=q1/q; ff=r_2*(nn-n-1)/(1-r_2)*n); r=sqrt(n*f/(nn-n-1+n*f)第一题：一：运行图形和结果：1969、1995和2000年的人口为：左图的计算为三次多项式拟合的函数计
23、算而得的结果：单位（亿）1969年人口为：g1969_3 =8.10404781077523 1969年人口为：g1995_3 =11.0391780295176 1969年人口为：g2000_3 =10.7846676285844 偏差平方和为：sumerror3 =0.472697009288482 右图的计算为二次多项式拟合的函数计算而得的结果：单位（亿）1969年人口为：g1969_2 =8.051768373432421969年人口为：g1995_2 =12.744093806231969年人口为：g2000_2 =13.8025210025012偏差平方和为：sumerror2 =
24、 0.709032276353744二：源代码（注：把t的值存放在t.txt中，p的值存放在p.txt中，并且把它们和M文件放到同一个文件夹里）clearx=load(x.txt);y=load(y.txt);n=length(x);c(1:n)=1; A=c,x,x.2,x.3;Q,R=qr(A,0);a=R(Q*y);a=fliplr(a); A1=c,x,x.2;Q,R=qr(A1,0);a1=R(Q*y);a1=fliplr(a1); t = 1949:2000; % 三次多项式拟合g = polyval(a,t);plot(x,y,o,t,g,m)hold ony1=a(1)*x.3
25、+a(2)*x.2+a(3)*x+a(4);g1969_3 = polyval(a,1969)g1995_3= polyval(a,1995)g2000_3 = polyval(a,2000)error3=y-y1;sumerror3=error3*error3 g = polyval(a1,t); %二次多项式拟合plot(x,y,o,t,g,r)hold ony2=a1(1)*x.2+a1(2)*x+a1(3);g1969_2 = polyval(a1,1969)g1995_2= polyval(a1,1995)g2000_2 = polyval(a1,2000)error2=y-y2;s
26、umerror2=error2*error2第二题：一：图形及结果：运算结果如下：以上五种拟合方法它们的均方误差如下：sum1 为指数函数拟合 y=a*exp(x)+b 的均方误差。sum2为三次多项式拟合的均方误差。sum3为二次多项式拟合的均方误差。sum4为y2=a*x2+b*x+c 拟合的均方误差。sum 为指数 y=a*exp(b*x) 拟合的均方误差。从均方误差可以看出来，三次多项式的拟合效果较好！二：源代码：（注：把t的值存放在t.txt中，p的值存放在p.txt中，并且把它们和M文件放到同一个文件夹里）cleart=load(t.txt);p=load(p.txt);%-指数函
27、数拟合 y=a*exp(x)+bm=length(t);c(1:m)=1;A1=c,exp(t);Q,R=qr(A1,0);a1=R(Q*p);x=-69.7:100;y=a1(1)+a1(2)*exp(x);axis(-80,100,-500,4000);plot(x,y,m);hold ony1=a1(1)+a1(2)*exp(t);error1=p-y1;sum1=sqrt(error1*error1)/m)%-三次多项式拟合A2=c,t,t.2,t.3;Q,R=qr(A2,0);a2=R(Q*p);a2=fliplr(a2);x=-69.7:100;y = polyval(a2,x);
28、axis(-80,100,-500,4000);plot(x,y,r)y1=a2(1)*t.3+a2(2)*t.2+a2(3)*t+a2(4);error2=p-y1;sum2=sqrt(error2*error2)/m) %-二次多项式拟合A3=c,t,t.2;Q,R=qr(A3,0);a3=R(Q*p);a3=fliplr(a3);x=-69.7:100;y = polyval(a3,x);axis(-80,100,-500,4000);plot(x,y,b)hold ony1=a3(1)*t.2+a3(2)*t+a3(3);error3=p-y1;sum3=sqrt(error3*err
29、or3)/m)%- y2=a*x2+b*x+c 拟合y=p.2;x=t;c(1:m)=1;A4=c,x,x.2;Q,R=qr(A4,0);a4=R(Q*y);x=-69.7:100;y=sqrt(a4(3)*x.2+a4(2)*x+a4(1);axis(-80,100,-500,4000);plot(t,p,o,x,y,black);y1=sqrt(a4(3)*t.2+a4(2)*t+a4(1);error4=y1-p;sum4=sqrt(error4*error4)/m)%-指数 y=a*exp(b*x) 拟合y=log(p);x=t;A=c,x;Q,R=qr(A,0);a=R(Q*y);x
30、=-69.7:100;a(1)=exp(a(1);y=a(1)*exp(a(2)*x);axis(-80,100,-500,4000);plot(t,p,o,x,y,g);hold ony1=a(1)*exp(a(2)*t);error=y1-p;sum=sqrt(error*error)/m)附加题第一题二次曲面拟合地表（结合自己专业知识）背景：一般地表是连续的曲面组成的，我们可以根据一些地面已知点（知道三维坐标）用曲面来拟合地表。步骤：首先：我们可以到实地去取样点，就是所说的去实地测量地表点的（X，Y，Z）坐标。在地表均匀地选取一些点进行三维坐标（X，Y，Z），然后我们可以根据所测量的点来
31、拟合地表曲面。其次：根据地表模型，一般用二次曲面来拟合地表，选择拟合公式。Z=f（X，Y）=a1 + a2*X + a3*Y + a4*XY + a5X2 + a6*Y2再次：运用数值分析中的曲线拟合知识。采样点：（x1,y1,z1）(x2,y2,z2)(xm,ym,zm) 基函数：（，一般）函数族：记号：，在这里： 采样点：我们已经实地测量出（X，Y，Z）值。 基函数：由 Z=f（X，Y）=a1 + a2*X + a3*Y + a4*XY + a5X2 + a6*Y2 我们6 个基函数：=1，=x =y=x*y =x2 =y2 我们可以求出 A ， 我们用Cholesky分解法求解，但直接计
32、算是极其数值不稳定的，不宜采用。最常用的方法是正交三角分解法（简称QR分解），它间接地对进行了Cholesky分解，其数值稳定性相当好。这里我们借用Matlab分解命令来求解。方法如下：得到，其中的列是规范正交的，即，是非奇异上三角矩阵。此时（间接做了Cholesky分解），法方程为这样只需求解上面上三角方程组，可用这样就求出了二次曲面函数的系数。最后：二次曲面函数绘图。调用MATLAB里的mesh（x,y,z）绘出曲面。绘出的图形如下：注：已知的（x,y,z）坐标是编的！没有进行实地测量！程序代码如下：clearx=load(x1.txt);y=load(y1.txt);z=load(z1.
33、txt);m=length(x);c(1:m)=1;%-用二次曲面拟合地表%-z=a(1)+a(2)*x+a(3)*y + a(4)*x.*y+a(5)*x.2+a(6)*y.2 A=c,x,y,(x.*y),x.2,y.2; Q,R=qr(A,0);a=R(Q*z); T = 1:40;S= 1:40;t,s=meshgrid(T,S);g=a(1)+a(2)*t+a(3)*s + a(4)*t.*s+a(5)*t.2+a(6)*s.2 ; mesh(t,s,g) 注：以上程序中 X1=-6 1 4 2 6 9 12 15 16 10 7 3 8 12 15 11 17 20 14 21 2
34、5 24 29 32 28 35 38 30 Y1=1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 Z1=25 27 26 29 35 28 27 15 29 38 46 50 26 20 32 16 29 35 40 56 50 48 43 40 35 31 28 20运行代码之前，把编好的x坐标放在x1.txt，y坐标放在y1.txt，z坐标放在z1.txt并把它们和M文件放在同一个文件夹下。然后就可以运行。第六章题目：求积分方程的数值解，并与精确解比较（精确解为）解：1、理论分析设离散化
35、为，下面求的近似值。对方程中的积分用某个求积公式：近似代替得：记得：这样转化为关于的线性方程组：其中是单位矩阵，2、实验结果及分析采用五等分和七等分的高斯柯特斯公式得到的结果如图1 图1 五、七等分的高斯柯特斯公式积分结果图五等分的高斯柯特斯公式运行结果和精确值分别为：0.99997600726849 1.22137345337173 1.49178890469185 1.82207508278336 2.22548753168652 2.71821660945298精确值为：1.00000000000000 1.22140275816017 1.49182469764127 1.822118
36、80039051 2.22554092849247 2.71828182845905七等分的高斯柯特斯公式运行结果和精确值分别为：0.99999987238056 1.15356484767778 1.33071202762260 1.53506281335132 1.77079472644729 2.04272680957445 2.35641814165885 2.71828148155343精确值为：1.00000000000000 1.15356499489511 1.33071219744735 1.535063009255211.77079495243515 2.04272707026614 2.35641844238366 2.71828182845905八等分的高斯柯特斯公式运行结果和精确值分别为：0.99999999923687 1.13314845220208 1.28402541570786 1.45499141350785 1.64872126944194 1.86824595600650 2.11700001499712 2.39887529213644 2.71828182638464精确值为：1.00000000000000 1.1331484
