计算物理学pptchp4

1、目录 实验数据的拟合计算物理 (Computational Physics)上节回顾实验数据的插值线性插值 拉格朗日线性插值(一次基函数、一次插值多项式) 牛顿线性插值(一阶差商、一次插值多项式)二次插值 拉格朗日二次插值(二次基函数、二次插值多项式) 牛顿二次插值(二阶差商、二次插值多项式) n 次插值n n 次拉格朗日插值（n次插值基函数，n次插值多项式） n 次牛顿插值 ( 1n各阶差商，n次插值多项式）逐次线性插值二元函数的拉格朗日多点插值计算物理第四章 实验数据的拟合计算物理 (Computational Physics)目录 实验数据的拟合计算物理 (Computational P

2、hysics)目录实验数据的拟合最小二乘法与一元线性拟合多元线性拟合非线性曲线拟合目录 实验数据的拟合计算物理 (Computational Physics)目录 实验数据的拟合 插值是这样一种手段：给出一组离散点，确定一个通过这些点的函数来近似原函数。在实际中，数据不可避免的会有误差，插值函数会将这些误差也包括在内。因此，我们需要一种新的手段来近似原函数：不要求过所有的点（可以消除误差影响）；尽可能表现数据的趋势，靠近这些点。 有时候，问题本身不要求构造的函数过所有的点。如：有5个景点，要修一条公路S使得S为直线，且到所有景点的距离和最小。实验数据的拟合计算物理 (Computational

3、 Physics)目录 实验数据的拟合拟合问题：通过已知数据估计物理量之间的函数关系。对于两个物理量 x 和 y，它们之间存在某种关系: y = f (x)。根据测量得到一组 x 和 y 的数据：拟合问题就是寻找函数 f (x) 的一个近似表达式 y(x), 但是 y(x)不一定经过上述测量值。xx1x2xnyy1y2yn实验数据的拟合计算物理 (Computational Physics)实验数据的拟合目录 实验数据的拟合 x=0:0.5:10; y=sin(x);分别用五次和六次多项式进行拟合。例：012345678910-1-0.500.511.5012345678910-1-0.500

4、.511.5计算物理 (Computational Physics)最小二乘法与一元线性拟合目录 实验数据的拟合金属电阻是温度的函数，设函数关系为R = R0(1+t)其中 R0 为 0 时的金属电阻, 为金属的电阻温度系数, 如何根据一组测量值求 R0 和 ?设一元拟合方程为：y(x) = A0 + A1x测量值例：xx1x2xnyy1y2yn计算物理 (Computational Physics)最小二乘法与一元线性拟合目录 实验数据的拟合将测量值代入拟合方程 y(x) = A0 + A1x , 就可以写出 n 个线性方程式：由此方程组可以求得A0和A1的很多组解。哪一组最接近真实值？哪一

5、组最接近真实值？计算物理 (Computational Physics)最小二乘法与一元线性拟合目录 实验数据的拟合已知一个函数的数值表已知一个函数的数值表，求一个求一个简单易算简单易算的近似函数的近似函数 y(x) f (x)。但是但是(1) n 通常很大通常很大；(2) yi 本身是测量值，不准确，即本身是测量值，不准确，即 yi f (xi) 。 这时没必要使这时没必要使 y(xi) = yi , 而只要而只要 yi y(xi) 总体上总体上尽可能小。尽可能小。常见做法：常见做法：使 最小最小二乘法最小二乘法xx1x2xnyy1y2yn计算物理 (Computational Physic

6、s)最小二乘法与一元线性拟合目录 实验数据的拟合令：令：求使平方和求使平方和 达到最小的达到最小的A0和和A1 即即:要求：要求：解得：解得：其中其中xx1x2xnyy1y2yn拟合方程 y(x) = A0 + A1x计算物理 (Computational Physics)目录 实验数据的拟合最小二乘法线性拟合的算法流程v计算v计算 x=0:10:200; y=4.38 4.56 4.70 4.86 5.08 5.24 5.40 5.58 5.74 5.96 6.06 6.26 6.44 6.58 6.74 6.94 7.12 7.28 7.42 7.60 7.78; A1=(mean(x.*

7、y)-mean(x)*mean(y)/(mean(x.2)-mean(x)2)A1 = 0.0170 A0=mean(y)-A1*mean(x)A0 = 4.3824 t=0:200; R=A0+A1*t; plot(t,R,r); hold on plot(x,y,.,Markersize,15)02040608010012014016018020044.555.566.577.58计算物理 (Computational Physics)最小二乘法与一元线性拟合目录 实验数据的拟合注意：（1）数据拟合必须首先确定拟合方程的形式。如果已知物理量之间存在某种形式的函数关系，只是由于偶然误差而造成

8、数据离散，这时拟合方程与原函数关系具有相同的形式。（2）如果两物理量之间关系式的形式尚不清楚，这时需要先用观测数据在以物理量为坐标的平面上画出散点图，再按照散点分布的大致规律，粗略地分析判定拟合方程的形式，然后才能使用前面讨论的拟合方法。 计算物理 (Computational Physics)最小二乘法与一元线性拟合目录 实验数据的拟合有些非线性函数关系可经过适当变换后用线性拟合求解。如半导体的电阻与温度的关系为:其中：k 是波耳兹曼常量, Eg 是禁带宽度, C 在温度范围不大时可视为常量。可见电阻 R 和温度T 是非线性关系，但对上式两边取对数后：InR与1/T呈线性关系，转化为线形拟合

9、问题。只须令：计算物理 (Computational Physics)最小二乘法与一元线性拟合目录 实验数据的拟合已知电流和电压的函数关系式为 I = aebV, 试用最小二乘法拟合此函数。Vi 1 2 3 4 5 6 7 8Ii15.3 20.5 27.4 36.6 49.1 65.6 87.8 117.60123456789020406080100120VI例： 某一物理实验中测得电压和电流的一组数据：计算物理 (Computational Physics)最小二乘法与一元线性拟合目录 实验数据的拟合将 I = aebV 化为线形模型 In I = In a + bV v=1:8; i=1

10、5.3 20.5 27.4 36.6 49.1 65.6 87.8 117.6; x=v; y=log(i); b=(mean(x.*y)-mean(x)*mean(y)/(mean(x.2)-mean(x)2)b = 0.2912 a=exp(mean(y)-b*mean(x)a = 11.4371 vv=0:0.1:9; ii=a*exp(b*vv); plot(vv,ii,r) plot(v,i,.,Markersize,15) hold on plot(vv,ii,r)0123456789020406080100120140160计算物理 (Computational Physics)

11、多元线性拟合目录 实验数据的拟合设物理量 y 随 k 个自变量 x1, x2, , xk 的变化而变化, 即 y = f (x1, x2, , xk)。 已知 n 组实验数据, x1i, x2i, , xki, yi (i=1,2,n), 试用多元线性函数拟合y = f (x1, x2, , xk)。设：最小二乘原则计算物理 (Computational Physics)多元线性拟合目录 实验数据的拟合令：则有即：最小二乘问题写成矩阵形式：计算物理 (Computational Physics)多元线性拟合目录 实验数据的拟合求最小值，则：计算物理 (Computational Physics

12、)多元线性拟合目录 实验数据的拟合对于一元线性拟合: i=15.3 20.5 27.4 36.6 49.1 65.6 87.8 117.6; Y=log(i)Y = 2.7279 3.0204 3.3105 3.6000 3.8939 4.1836 4.4751 4.7673 v=1:8; X=ones(8,1) v%空格连接两个列向量X = 1 1 1 2 1 3 1 4 1 5 1 6 1 7 1 8 Xans = 1 1 1 1 1 1 1 1 1 2 3 4 5 6 7 8 A=(X*X)(X*Y)A = 2.4369 0.2912 a=exp(A(1)a = 11.4371 b=A(

13、2)b = 0.2912计算物理 (Computational Physics)多元线性拟合目录 实验数据的拟合设拟合函数为： y(x1, x2) = A0 + A1 x1 + A2 x2例：设 y 是 x1 和 x2 的线性函数，已知下列一组测量值，试拟合二元一次线性函数：x10.1934 0.6822 0.3028 0.5417 0.1509 0.6979 0.3784 0.8600 0.8537 0.5936x20.4966 0.8998 0.8216 0.6449 0.8180 0.6602 0.3420 0.2897 0.3412 0.5341y2.2455 3.5731 2.950

14、5 2.9075 2.7214 3.1143 2.0166 2.3899 2.6112 2.7092计算物理 (Computational Physics)多元线性拟合目录 实验数据的拟合 x1=0.1934 0.6822 0.3028 0.5417 0.1509 0.6979 0.3784 0.8600 0.8537 0.5936; x2=0.4966 0.8998 0.8216 0.6449 0.8180 0.6602 0.3420 0.2897 0.3412 0.5341; y=2.2455 3.5731 2.9505 2.9075 2.7214 3.1143 2.0166 2.3899

15、 2.6112 2.7092; X=ones(10,1) x1 x2; Y=y; A=(X*X)(X*Y)A = 0.8930 1.1280 2.1173 X1,X2=meshgrid(0:0.05:1); Y0=A(1)+A(2)*X1+A(3)*X2; surf(X1,X2,Y0) alpha(0.6) hold on plot3(x1,x2,y,.)计算物理 (Computational Physics)目录 实验数据的拟合三维图形矩阵网格生成函数 meshgrid(), 使用格式：X,Y = meshgrid(x,y)，分别产生以向量x为行，向量y为列的两个大小相同的矩阵。X,Y =

16、meshgrid(x),等同于meshgrid(x,x)。 X,Y=meshgrid(0:4)X = 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4Y = 0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 X,Y=meshgrid(0:3,1:0.1:1.4)X = 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3Y = 1.0000 1.0000 1.0000 1.0000 1.1000 1.1000 1.1000 1.1000 1.2000 1.2000 1.

17、2000 1.2000 1.3000 1.3000 1.3000 1.3000 1.4000 1.4000 1.4000 1.4000计算物理 (Computational Physics)目录 实验数据的拟合曲面绘制函数plot3(), mesh(), meshc(), meshz(), surf(), surfc(), surfl()plot3() 绘制三维曲线图：z = x2 + y2 X,Y=meshgrid(-2*pi:0.2:2*pi); Z=X.2+Y.2; plot3(X,Y,Z) t=0:0.1:8*pi; plot3(sin(t),cos(t),t)-10-50510-10

18、-50510020406080-1-0.500.51-1-0.500.51051015202530计算物理 (Computational Physics)目录 实验数据的拟合曲面绘制函数plot3(), mesh(), meshc(), meshz(), surf(), surfc(), surfl()mesh() 绘制三维网格图： X,Y=meshgrid(-3*pi:0.4:3*pi); R=sqrt(X.2+Y.2); Z=sin(R)./R; mesh(X,Y,Z)-10-50510-10-50510-0.500.51-10-50510-10-50510-0.500.51hidden

19、off：表示显示遮住部分hidden on：表示不显示遮住部分计算物理 (Computational Physics)目录 实验数据的拟合曲面绘制函数plot3(), mesh(), meshc(), meshz(), surf(), surfc(), surfl()meshc() 生成具有基本等高线的网格图meshz() 生成带有基准平面的网格图 meshc(X,Y,Z) meshz(X,Y,Z)-10-50510-10-50510-0.500.51-10-50510-10-50510-0.500.51计算物理 (Computational Physics)目录 实验数据的拟合曲面绘制函数p

20、lot3(), mesh(), meshc(), meshz(), surf(), surfc(), surfl()surf() 生成三维表面图surfc() 生成具有等高线的表面图surfl() 生成具有光照效果的表面图 surf(X,Y,Z) surfc(X,Y,Z) surfl(X,Y,Z)-10-50510-10-50510-0.500.51-10-50510-10-50510-0.500.51-10-50510-10-50510-0.500.51计算物理 (Computational Physics)目录 实验数据的拟合shading flat：去掉线条 surf(X,Y,Z) sh

21、ading flat-10-50510-10-50510-0.500.51shading interp: 去掉线条,平滑颜色 surf(X,Y,Z) shading interpalpha() 设置透明度, 参数在0 1上。 surf(X,Y,Z) shading interp alpha(0.6)计算物理 (Computational Physics)目录 实验数据的拟合sphere(): 生成单位球体X,Y,Z = sphere(n); X,Y,Z,为(n+1)(n+1)维矩阵 X,Y,Z = sphere(30); surf(X,Y,Z) axis equal mesh(2*X,2*Y,

22、2*Z) axis equal hidden off-1-0.500.51-0.500.5-1-0.500.51-2-1012-101-2-1012 mesh(X,Y,2*Z) axis equal-101-0.500.5-2-1.5-1-0.500.511.52计算物理 (Computational Physics)目录 实验数据的拟合X,Y,Z = cylinder(r,n) 画高度为1的圆柱体r 为半径, n为图形的拼片个数 X,Y,Z = cylinder(0.5,30); surf(X,Y,Z) axis equal-0.500.5-0.4-0.200.20.400.20.40.60

23、.81 t = 0:pi/10:2*pi; X,Y,Z = cylinder(2+cos(t),30); surf(X,Y,Z)-4-2024-4-202400.20.40.60.81计算物理 (Computational Physics)目录 实验数据的拟合 X,Y,Z = cylinder(1+cos(t),30); surf(X,Y,Z)-2-1012-2-101200.20.40.60.81 t = 0:pi/10:4*pi; X,Y,Z = cylinder(1+cos(t),30); surf(X,Y,4*pi*Z) axis equal-202-101024681012计算物理

24、(Computational Physics)非线性曲线拟合目录 实验数据的拟合用 k (k n) 次多项式拟合, 拟合函数为：y(x) = A0 + A1 x + A2 x2 + + Ak xk最小二乘原则：测量值：xx1x2xnyy1y2yn计算物理 (Computational Physics)非线性曲线拟合目录 实验数据的拟合矩阵形式：解得：令：计算物理 (Computational Physics)目录 实验数据的拟合function f=ls(x,y,k)% 最小二乘曲线拟合,生成k次多项式sx=size(x); %size( ) 显示矩阵大小if sx(1)=1 x=x;ends

25、y=size(y);if sy(1)=1 y=y;endn=max(sx);if length(y)=n error(Vectors must be the same lengths);elseif (n-1) x=0:10; y=cos(x); ls(x,y,5);012345678910-1.5-1-0.500.511.5计算物理 (Computational Physics)目录 实验数据的拟合size( ) 显示矩阵大小显示矩阵大小 A=1 2 3;4 5 6; size(A)ans = 2 3error() 程序到此结束程序到此结束,显示引号中内容显示引号中内容 x=0:10; y=

26、1 2 3; ls(x,y,5);? Error using = lsVectors must be the same lengthsrepmat(), 由向量生成矩阵由向量生成矩阵 x=1 2 3; repmat(x,2,1)ans = 1 2 3 1 2 3cumprod( ), 生成连积矩阵生成连积矩阵B = cumprod(A, dim)此调用格式返回A中由标量dim所指定的维数的累加积: A=1 2 3;4 5 6A = 1 2 3 4 5 6 cumprod(A)%沿着第一维（列）的累加积ans = 1 2 3 4 10 18 cumprod(A,2)%沿着第二维（行）的累加积ans = 1 2 6 4 20 120计算物理 (Computational Physics)非线性曲线拟合目录 实验数据的拟合Polyfit 多项式拟合命令格式： p=polyfit(x,y,n)其中：p为返回的多项式向量，从最高次幂到最低次幂的系数x为准备拟合的自变量，y为应变量n为拟合的次数计算物理 (Computational Physics)非线性曲线拟合目录 实验数据的拟合例：cos(x)的曲线拟合x=0:10;y=cos(x);a=polyfit(x,y,5) ;x1=0:0.05:10;y1=a(6)+a(5)*x1+