第四章 数值计算.ppt

1、第四章 数值计算,4.1 解线性方程组 4.2 曲线拟合 4.3 数据差值,4.1 解线性方程,MATLAB解线性方程一般不用先求逆阵，再做矩阵乘法的方法求解。而是采用类似数量除法的方法来描述一般的线性方程组的解。 系数矩阵左乘未知数 Ax=B x=BA 系数矩阵右乘未知数 xA=B x=B/A 系数矩阵为m行n列矩阵 (1) m=n 为方阵系统，可以求出一个精确的解。 (2) mn 过限定系统，可以求出一个最小方阵的解。 (3) mn 限定不足系统，可以求出一个至多含有m个非零元素的基解。,2.6.1 方阵系统 Eg: A= 2 3 4 Eg: A= 2 3 4 5 4 1 5 4 1 1

2、3 2 1 3 2 U= 1 B= 4 7 5 2 1 4 7 3 5 2 6 Ax=U x=A-1U=AU Ax=B x=A-1B=AB x= -0.8148 x= -1.4444 1.2222 -0.4074 1.6667 2.0000 -1.0000 2.3333 -0.5926 0.2222 1.8889 -0.2963 A*x A*x ans = 1 ans= 4 7 5 2 1 4 7 3 5 2 6,2.6.2 过限定系统 各种曲线拟合的过程中经常遇到过限定系统的线性方程组，MATLAB将两组数据保存在列矢量t与列矢量y中。 Eg: t=0 0.3 0.8 1.1 1.6 2.3

3、 y=0.82 0.72 0.63 0.60 0.53 0.50 认为这组数据大概在如下一个衰减的指数函数曲线上y(t)=C1+C2e-t y是两个矢量的线性组合，一个是常数C1，另一个是e-t，将数据代入函数y(t)=C1+C2e-t得到六个含有两个未知数的等式，将这些数据表示成一个六行两列的矩阵E，把未知参数看作一个矢量c的两个元素方程组为Ec=y。,解： t=0 0.3 0.8 1.1 1.6 2.3； y=0.82 0.72 0.63 0.60 0.53 0.50； E=ones(size(t)exp(-t) E =1.0000 1.0000 1.0000 0.7408 1.0000

4、0.4493 1.0000 0.3329 1.0000 0.2019 1.0000 0.1003 c=Ey c = 0.4760 0.3413 E*c ans= 0.8173 0.7288 0.6293 0.5896 0.5449 0.5102,4.2 曲线拟合,一回归法曲线拟合 曲线拟合需要从一些离散的数据中推导出两者之间的数学解析关系，而数据插值是通过原始数据计算一些新的离散数据点。曲线拟合的结果一般为一个或者多个数学解析关系，利用这些解析关系能够对数据进行一定的推断，准确的曲线拟合结果可以用来进一步评估、验证实测的数据。 利用MATLAB进行曲线拟合主要有两种方法：回归法拟合和多项式拟合

5、。,t=0 0.3 0.8 1.1 1.6 2.3 y=0.5 0.82 1.14 1.25 1.35 1.40 采用多项式 y=a0+a1t+a2t2 解线性方程 Y=XA A=XY x=ones(size(t) t t.2 x=1 0 0 1 0.3 0.09 1 0.8 0.64 1 1.1 1.21 1 1.6 2.56 1 2.3 5.29,a=xY a= 0.5318 这个二次多项式函数为 y=0.5318+0.9191t-0.2387t2 0.9191 - 0.2387 T=(0:0.1:2.5) ; Y=ones(size(T) T T.2*a plot(T,Y,-,t,y,o

6、); grid on,2.采用一个带有线性参数的指数函数 y=a0+a1e-t+a2te-t x=ones(size(t) exp(-t) t.*exp(-t); a=xY a= 1.3974 -0.8988 0.4097 因此数据的拟合函数为 y=1.3974-0.8988e-t+0.4097te-t; T=(0:0.1:2.5) ; Y=ones(size(t) exp(-T) T.*exp(-T)*a; plot(T,Y,-,t,y,o); grid on 例： t = 0 .3 .8 1.1 1.6 2.3; y = 0.5 0.82 1.14 1.25 1.35 1.40; plot

7、(t,y,r*) grid on, X = ones(size(t) t t.2 X = 1.0000 0 0 1.0000 0.3000 0.0900 1.0000 0.8000 0.6400 1.0000 1.1000 1.2100 1.0000 1.6000 2.5600 1.0000 2.3000 5.2900 A = Xy A = 0.5318 0.9191 -0.2387 这样得到的多项式应该为 y=0.5318+0.9191t-0.2387t2 为了验证结果，可以进一步进行计算： T = (0:0.1:2.5); Y = ones(size(T) T T.2*A; plot(T,

8、Y,-,t,y,o), grid on legend(Fitting,Origin),除了多项式以外，还可以猜测原始数据的多项式是由指数函数组成的： 于是，为了求得系数，可以在MATLAB的命令行中键入下面的指令 X = ones(size(t) exp(-t) t.*exp(-t) X = 1.0000 1.0000 0 1.0000 0.7408 0.2222 1.0000 0.4493 0.3595 1.0000 0.3329 0.3662 1.0000 0.2019 0.3230 1.0000 0.1003 0.2306 左除： A = Xy A = 1.3974 -0.8988 0.

9、4097 这样得到的多项式应该为 为了验证这次得到的结果，可以进一步计算： T = (0:0.1:2.5); Y = ones(size(T) exp(-T) T.*exp(-T)*A; plot(T,Y,-,t,y,o), grid on legend(Fitting, Origin),二、多项式曲线拟合 根据一组已知的自变量和函数的值，应用最小二乘法求出多项式的拟合曲线 polyfit p=polyfit(x,y,n) x 已知自变量的取值 y 相应的函数值 n 为多项式的次数 Eg: x=1 2 3 4 5 y=6.5 43.1 128 290.7 498.4 求一个三次多项式与这组数据

10、拟合 p=polyfit(x,y,3) x2=1:0.1:5 y2=polyval(p,x2) plot(x,y,o,x2,y2) grid on,001%CURVE_FIT 多项式拟合计算示例 002% 准备原始数据 003x = 0:.1:10; 004y = sin(x)+cos(2*x); 005% 5次多项式拟合 006k5 = polyfit(x,y,5); 007y5 = polyval(k5, x); 008% 11次多项式拟合 009k11 = polyfit(x,y,11); 010y11 = polyval(k11, x); 011% 绘制数据曲线 012plot(x,

11、y, g, x, y5, r, x, y11, b) 013% 标注 014title(Curve Fitting, Fontsize,14) 015legend(Original Curve, 5th order, 11th order,4 ) 016set(findobj(Type, line), LineWidth, 2),4.3 数据插值,插值就是在已知数之间寻找估计值。,1、多项式插值 格式：yi=interp1(x,y,xi,method) xi用于插值的点 method指定插值的方法，* method可以为插值计算指定相应的算法，为字符串类型，其取值可以为nearest、line

12、ar、spline、cubic、pchip、v5cubic。 插值方法: （1）邻近插值：在已知数据最邻近的点设置插值。 （2）线性插值(默认)：在已知的成对数据点与离散函数之间进行拟合，返回xi指定的相应的数值。 （3）三次样条插值：采用一个三次函数对已知数据进行拟合，返回xi的函数值。,Eg: t表示19001990年人口普查年，每10年普查一次，p为当年人口数量 t=1900:10:1990 p=75.995 91.972 105.711 123.203 131.669 150.697 179.323 203.212 226.505 249.633 interp1(t,p,1975) a

13、ns=214.8585 x=1900:1:2000 y=interp1(t,p,x,spline); plot(t,p,o,x,y),001%INTERP_EX1 一维插值计算示例 002% 准备数据 003x = 0:10; 004y = cos(x); 005% 插值点 006xi = 0:0.2:10; 007% 进行插值运算 008yin = interp1(x,y,xi,nearest); 009yic = interp1(x,y,xi,cubic); 010% 绘制结果 011plot(x,y,o,xi,yin,*,xi,yic) 012legend(origin,nearest,

14、cubic) 013title(一维插值计算示例),001function compare_interp( ) 002%COMPARE_INTERP 不同插值运算的比较 003 004% 原始数据 005x,y = meshgrid(?3:1:3); 006z = peaks(x,y); 007figure(1); clf 008surfc(x,y,z); 009title(原始数据) 010% 进行插值运算 011xi, yi = meshgrid(?3:0.25:3); 012zi1 = interp2(x,y,z,xi,yi,nearest); 013zi2 = interp2(x,y,

15、z,xi,yi,linear); 014zi3 = interp2(x,y,z,xi,yi,cubic); 015zi4 = interp2(x,y,z,xi,yi,spline); 016% 通过可视化结果比较 017figure(2) 018subplot(2,2,1);surf(xi,yi,zi1);,019title(二维插值 - nearest) 020subplot(2,2,2);surf(xi,yi,zi2); 021title(二维插值 - linear) 022subplot(2,2,3);surf(xi,yi,zi3) 023title(二维插值 - cubic) 024subplot(2,2,4);surf(xi,yi,zi4) 025title(二维插值 - spline) 026% 可视化结果 027figure(3) 028