Matlab详细教程

HUASHI3483

7.3.2 多项式回归  

------------------------------------------------------------------------------
--  

当使用 polyfit 时的 的回归法，我们统称这些为多项式回归。不过值得注意的是，不见
得用高阶的多项式就能合理的代表原数据，此点可以由上图来说明。图中所做的曲线契合
分别从二阶到九阶的多项式。从图中可以看出越高阶的多项式所形成的方程式的振荡程度
越剧烈（七阶以上的皆有此现象），另外五阶以上的多项式都会通过所有的原始数据点。
至于选择几阶的多项式是合理的曲线契合函数，就见仁见智了。   



------------------------------------------------------------------------------
--  
   

HUASHI3483

7.3.3 多项式契合及函数计算  

------------------------------------------------------------------------------
--  

和 polyfit有关的另一个函数是用来做多项式函数计算 polyval，由 polyfit 算出多项
式的各个系数 ()后，即可以 polyval 计算多项式的函数值。它的语法为
polyval(coef,x)，其中 coef 是多项式的各个系数所构成的阵列，x则是要计算多项式值
的x阵列。举一例说明   

＞> x=[0 1 2 3 4 5];   

＞> y=[0 20 60 68 77 110];   

＞> coef=polyfit(x,y,1); % 计算线性回归的各项系数   

＞> ybest=polyval(coef,x); % 直接以polyval计算多项式的数值   


这个函数 polyval真是大大的好用，如果没有它的话，那我们要计算多项式的函数值还要
费一番功夫。假设已有五阶的多项式回归，要计算对应的多项式值，则须以下的步骤   

＞> coef=polyfit(x,y,5);   

＞> a0=coef(1);   

＞> a1=coef(2);   

＞> a2=coef(3);   

＞> a3=coef(4);   

＞> a4=coef(5);   

＞> a5=coef(6);   

＞> f=a0 + a1*x + a2*x.^2 + a3*x.^3 + a4*x.^4 + a5*x.^5;  


我们接著示范从二阶到九阶的多项式回归的程式   

＞> x=[0 1 2 3 4 5];   

＞> y=[0 20 60 68 77 110];   

＞> newx=0:0.05:5; % 以新阵列形成更小增量以利回归计算及绘图   

＞> for n=2:9   

＞> f(:,n)=polyval(polyfit(x,y,n),newx)';   

＞> plot(newx,f(:,n),x,y,'o')   

＞> title(['Poly. regression, deg=',int2str(n)])   

＞> xlabel('Time'), ylabel('Temp'), grid   

＞> pause % 每次要暂停看清楚图再执行下一步骤   

＞> end  


在上述的 title指令中我们示范了如何将一变数输入（n代表多项式的阶数），是利用
int2str这个指令，它是用来将一整数(integer) 转换成为一个字串 (string)，因为在
title中只能以字串出现。此外在 title指令中尚须以 [ ] 将所有叙述包括在内。类似的
指令尚有 num2str 它是用来将实数转换成一字串。有关这几个新介绍的详细说明，请参
考 title, int2str, num2str 的线上说明。   



-----------------------

HUASHI3483

8.1 多项式的根  

------------------------------------------------------------------------------
--  

一个多项式视其阶数而定，它的根可以有一个到数个，可能为实数也可能是复数。要求一
高阶多项式的根往往须借助数值方法，所幸MATLAB已将这些数值方法写成一函数roots(p)
，我们只要输入多项式的各阶系数（以 p 代表）即可求解到对应的根。   

＞> p=[1 3 2];   

＞> r=roots(p)   

r =   

-2   

-1   

＞> p=[1 -12 0 25 116]; % 注意二阶项系数为零须要输入，否则多项式的阶数就不对   

＞> r=roots(p) % 有实数根及复数根   

r =   

11.7473   

2.7028   

-1.2251 + 1.4672i   

-1.2251 - 1.4672i   


与 roots 相关的函数尚有 poly, real，这二个函数的用途是要验算求解的根展开能求得
原多项式。例如有一个二次方程式的根为2, 1，则以下式计算原多项式   

   

poly 函数就是在求出多项式的各阶系数，其语法为 poly(r)，其中 r 是代表根的阵列。
而 real 则是用来去除因计算时产生的假虚部系数，为何会有此种情形请参考以下的例子
。   

＞> r=[-2 1];   

＞> pp=poly(r) % pp=(x+2)(x-1)=x^2+3x+2   

pp =   

1 3 2   

＞> p=[1 -4 6 -4];   

＞> r=roots(p)   

r =   

2.0000 1.0000 + 1.0000i 1.0000 - 1.0000i   

＞> pp=poly(r) % 这个多项式的系数与原多项式 p 相同   

pp =   

1 -4 6 -4   

＞> pp=[1 7 12 9]; % 再看另一个多项式   

＞> r=roots(pp)   

r =   

-4.9395   

-1.0303 + 0.8721i   

-1.0303 - 0.8721i   

＞> pp=poly(r) % 注意因计算的误差会有假虚部产生   

pp =   

1.0000 7.0000 12.0000 9.0000 + 0.0000i   

＞> pp=real(pp) % 可以real将假虚部去除，将原多项式还原   

pp =   

1.0000 7.0000 12.0000 9.0000  



------------------------------------------------------------------------------
--  
   

HUASHI3483

8.2 非线性方程式的实根   

------------------------------------------------------------------------------
--  

如果求根的方程式不为多项式的形态 就不能用 roots 函数。而这类的方策’h半是非线
性方程式， 其函数形态变化很大。对于解这类方程式的根，可以用 fzero函数，它其实
是用来找一函数 f(x) 的 x 值代入时，会使该函数值为零 (f(x)=0)；而这也就是根的特
性，因此我们可以用 fzero求根。   


要求任一方程式的根有三步骤：   

先定义方程式。要注意必须将方程式安排成 f(x)=0 的形态，例如一方程式为sin(x)=3，
则该方程式应表示为 f(x)=sin(x)-3。可以 m-file 定义方程式。   
代入适当范围的 x, y(x) 值，将该函数的分布图画出，藉以了解该方程式的「长相」。  

由图中决定y(x)在何处附近(x0)与 x 轴相交，以fzero的语法fzero('function',x0) 即
可求出在 x0附近的根，其中 function 是先前已定义的函数名称。如果从函数分布图看
出根不只一个，则须再代入另一个在根附近的 x0，再求出下一个根。   
以下分别介绍几数个方程式，来说明如何求解它们的根。  


例一、方程式为   

sin(x)=0   

我们知道上式的根有 ，求根方式如下：   

＞> r=fzero('sin',3) % 因为sin(x)是内建函数，其名称为sin，因此无须定义它   

r = % 选择 x=3 附近求根   

3.1416   

＞> r=fzero('sin',6) % 选择 x=6 附近求根   

r =   

6.2832  


例二、方程式为先前提到的 MATLAB 内建函数 humps，我们不须要知道这个方程式的形态
为何，不过我们可以将它划出来，再找出根的位置。求根方式如下：   

＞> x=linspace(-2,3);   

＞> y=humps(x);   

＞> plot(x,y), grid % 由图中可看出在0和1附近有二个根   

＞> r=fzero('humps',1.2)   

r =   

1.2995  


例三、方程式为   

   

这个方程式其实是个多项式，我们说明除了用 roots 函数找出它的根外，也可以用这节
介绍的方法求根，注意二者的解法及结果有所不同。求根方式如下：   

% m-function, f_1.m   

function y=f_1(x) % 定义 f_1.m 函数   

y=x.^3-2*x-5;  


＞> x=linspace(-2,3);   

＞> y=f_1(x);   

＞> plot(x,y), grid % 由图中可看出在2和-1附近有二个根   

＞> r=fzero('f_1',2); % 决定在2附近的根   

r =   

2.0946   

＞> p=[1 0 -2 -5]   

＞> r=roots(p) % 以求解多项式根方式验证   

r =   

2.0946   

-1.0473 + 1.1359i   

-1.0473 - 1.1359i   


例四、方程式为   

   

求根方式如下：   

% m-function, f_2.m   

function y=f_2(x) % 定义 f_2.m 函数   

y=x.^2.*sin(x)+cos(x);   


＞> x=linspace(-3,3);   

＞> y=f_2(x);   

＞> plot(x,y), grid % 由图中可看出在-1和3附近有二个根   

＞> r=fzero('f_2',-1); % 决定在-1附近的根   

r =   

-0.8952   

＞> r=fzero('f_2',3); % 决定在3附近的根   

r =   

3.0333  


例五、方程式为   

   

求根方式如下：   

% m-function, f_3.m   

function y=f_3(x) % 定义 f_3.m 函数   

y=2*exp(-x).*sin(2*pi*x)-0.5;   


＞> x=0:0.1:2; y=f_3(x);   

＞> plot(x,y), grid % 由图中可看出在0,0.5和1附近有三个根   

＞> r=fzero('f_3',0) % 决定在0附近的根   

r =   

0.0420   

＞> r=fzero('f_3',0.5) % 决定在0.5附近的根   

r =   

0.4368   

＞> r=fzero('f_3',1) % 决定在1附近的根   

r =   

1.1435   


------------------------------------------------------------------------------
--  

     

HUASHI3483

9.1.1 梯形法  

------------------------------------------------------------------------------
--  

MATLAB提供最简单的积分函数是梯形法trapz，我们先说明梯形法语法trapz(x,y)，其中
x,y分别代表数目相同的阵列或矩阵，而y与x的关系可以由是一函数型态（如y=sin(x)）
或是不以函数描述的离散型态（像第八章介绍的x与y皆为离散点）。   


我们看一简单积分式   

   

以下为 MATLAB 的程式   

＞> x=0:pi/100:pi;   

＞> y=sin(x);   

＞> k=trapz(x,y)   

k =   

1.9998  



------------------------------------------------------------------------------
--  

HUASHI3483

9.1.2 二次函数法  

------------------------------------------------------------------------------
--  

MATLAB 另外提供二种积分函数，它们分别是辛普森法 quad 和牛顿-康兹法 quad8。三种
方法的精确度由低而高，分别为 trapz, quad, quad8。  


由于这二种方法依据的积分法不同于梯形法，因此它们的语法就和 trapz 不同；其语法
为 quad('function',a,b) （quad8语法相同），其中function是一已定义函数的名称（
如sin, cos, sqrt, log 等），而 a, b是积分的下限和上限。和 trapz比较，quad,  
quad8不同之处在于这二者类似解析式的积分式，只须设定上下限及定义要积分的函数；
而 trapz则是针对离散点型态的数据做积分。   


我们看一简单积分式   

   

以下为 MATLAB 的程式   

＞> a=0; b=0.5;   

＞> kq=quad('sqrt',a,b)   

kq =   

0.2357   

＞> kq8=quad8('sqrt',a,b)   

kq8 =   

0.2357  


再来看一个较复杂的积分式   

＞> x=-1:0.17:2;   

＞> y=humps(x);   

＞> area=trapz(x,y)   

area =   

25.9174   

＞> x=-1:0.07:2;   

＞> y=humps(x);   

＞> area=trapz(x,y)   

area =   

26.6243   

＞> area=quad('hump',-1,2)   

area =   

26.3450   

＞> area=quad8('hump',-1,2)   

area =   

26.3450  



------------------------------------------------------------------------------
--  
   

HUASHI3483

9.3.1 差分表示法  

------------------------------------------------------------------------------
--  

我们知道一微分项的计算，可以在二相邻点 x+h 和 x 间函数取下列极限求得   

   

若将原连续的空间以多个离散点取代，即是。则上述的极限以离散点的方式计算，即是以
下的差分式 (difference equation)   

   

其中，而上式被称为前向差分，因为是以为参考点，另一点在它之前。上式的几何意义是
以二点函数值计算在的斜率。事实上，除了上式可计算在的微分值外，也可以下列二式计
算   




而高阶微分项可以利用低阶微分项来计算，举例来说一个二阶微分式可以表示为   

   

所以对应的差分式有   





------------------------------------------------------------------------------
--  
   

HUASHI3483

9.3.2 差分函数  

------------------------------------------------------------------------------
--  

上述提及的后向差分式，在 MATLAB 有对应的diff 函数来计算二相临点的差值，它的语
法为 diff(x)，其中 x代表一组离散点 。假设有x, y(x)的数据为   

x=[1 3 5 7 9], y=[1 4 9 16 25]   

则 diff(x)=[2 2 2 2], diff(y)=[3 5 7 9]，注意二者皆以后向差分计算且数据点只剩  
4 个而不是5个。而的数值微分则为 dy=diff(y)./diff(x)。   


因此要计算下列多项式在 [-4, 5] 区间的微分   

   

可依以下方式求解   

＞> x=linspace(-4,5); % 产生100个x的离散点   

＞> p=[1 -3 -11 27 10 -24];   

＞> f=polyval(p,x);   

＞> plot(x,f) % 将多项式函数绘图   

＞> title('Fifth-deg. equation')   

＞> dfb=diff(f)./diff(x); % 注意要分别计算diff(f)和diff(x)   

＞> xd=x(2:length(x)); % 注意只有99个df值，而且是对应x2,x3,...,x100的点   

＞> plot(xd,dfb ) % 将多项式的微分项绘图   

＞> title('Derivative of fifth-deg. equation')  


假设我们要找出上述多项式的局部极值(local critical value)，可以利用局部极值所在
点的微分为零，所以在这个点的左侧及右侧的二点微分值一定是一正一负，因此这二点的
微分值相乘为负值，可针对此性质找出局部极值所在点。以原多项式的数据示范如下：  


＞> product=dfb(1:length(dfb)-1).*dfb(2:length(dfb)); % 注意二相临点值相乘的写
法   

＞> crit=xd(find(product<0)) % 注意用到find指令   


若上述多项式以中央差分方式计算其函数微分项，就不能以diff计算，而须自行计算如下
：   

＞> num=f(3:length(f))-f(1:length(f)-2); % 注意中央差分是 f(k+1)-f(k-1)   

＞> deno=x(3:length(f))-x(1:length(f)-2); % 注意中央差分是 x(k+1)-x(k-1)   

＞> df_c=num./deno;   

＞> xd=x(2:length(x)-1); % xd的点数只有98个   

＞> plot(xd,df_c)   

＞> title('Derivative of fifth-deg. polynomial')  


以下的例子是针对数据组为离散型态，要注意的是原数据所代表的函数分布并无明显的折
角，但是它的一阶微分经数值微分计算后的微分函数分布却有极大的曲折变化。   

＞> x=0:0.1:1;   

＞> y=[-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2];   

＞> plot(x,y,'o',x,y)   

＞> title('y(x) data plot')   

＞> ylabel('y(x)'), xlabel('x')   

＞> dy=diff(y)./diff(x);   

＞> xd=x(1:length(x)-1);   

＞> plot(xd,dy)   

＞> title('Approximate derivative using diff')   

＞> ylabel('dy/dx'), xlabel('x')  



------------------------------------------------------------------------------
--  

HUASHI3483

第十章 常微分方程式  

------------------------------------------------------------------------------
--  

工程问题挑战：载具性能  


在此章我们简单的说明什么是常微分方程式，以及如何来应用数值方法求解方程式。   


在这章中提到的工程问题挑战「载具性能」，所谓载具(vehicle) 是泛指汽车、飞机、太
空船。而这里谈的是飞机推进动力科技的一项进展，一种新式的无导管风扇(Unducted  
Fan, UDF) 的涡轮螺旋桨引擎(turboprop engine)。UDF的发展起源于1980年代中期，它
利用新的材料、创新的叶片外型和更高叶片转速，可使螺旋桨客机飞行时速和喷射客机 (
以涡轮风扇为动力)一样快，而又更省油。   



------------------------------------------------------------------------------
--  

第十章 常微分方程式   
10.1 微分方程式   
10.2 阮奇-库达方法   
10.3 范例问题：飞机发动机的加速性能分析   
10.4 高阶常微分方程式   

------------------------------------------------------------------------------
--  
   

HUASHI3483

10.1 微分方程式  

------------------------------------------------------------------------------
--  

一个一阶常微分方程式 (ordinary differential equation, ODE) 可以下式表示   

   

其中x 为独立变数而 y 是 x 的函数，我们就是要求解什么函数 y(x) 能满足上述 ODE。
以下有几个一阶 ODE 的例子：   

   

除了上述的已知ODE外，还须有起始条件y0=y(x0)才能解方程式，即是在x=x0时，y(x)=y0
。上述各个方程式的解析解 (analytical solution) 如下：   




以数值方法求解上述的 ODE 的问题，可以转换为在已知 y(a) 的函数值而要计算y(b)，
依据泰勒序数对 y(b) 做展开   

   

其中 b=a+h。一个一阶的泰勒序数近似式为   

   

而二阶的泰勒序数近似式为   




MATLAB 所依据解ODE 的数值方法就是利用像上述的二阶及更高阶的三、四、五阶泰勒序
数近似式来计算 f(b)。  



------------------------------------------------------------------------------
--