matlab数值运算课件

第四章数值计算功能4.1多项式计算4.2数值插值和曲线拟合4.3函数极值4.4非线性方程问题求解4.5常微分方程的数值求解4.6数值微分与积分4.7概率统计4.1多项式运算

1．多项式表示法2．多项式运算函数

多项式求根多项式求值多项式乘法和除法多项式的微积分1．多项式表示法

(1)直接法:多项式多项式系数用行向量表示，按降幂排列；缺某幂次项，其幂次项系数为零。P=[a0,a1,...an-1,an]符串形式,x多项式变量y=poly2str(P,'x')(2)字符串表示法：y=poly2sym(P,x)(3)符号多项式表示法：完整形式p=[1,3,0,4,5]y=poly2str(p,'x')y=x^4+3x^3+4x+5symsxy=poly2sym(p,x)y=x^4+3*x^3+4*x+5>>p=[1,3,0,4,5]p=13045>>y=poly2str(p,'x')y=x^4+3x^3+4x+5>>rp=roots(p)rp=-3.23460.5594+1.1980i0.5594-1.1980i-0.8843>>P=[1,8,3,4,-10];>>pp=poly2str(P,'X')pp=X^4+8X^3+3X^2+4X-10例2.多项式的运算函数多项式的值；根和微分；拟合曲线；部分分式polyvalpolyvalpolyvalmcovolution

（1）多项式求根（roots）（1）多项式的根=0

n次多项式n个根，或实根，或若干对共轭复根。

r=roots(P)P：多项式系数向量；r：n个根b(1),b(2),…,b(n)的向量。>>r=roots(P)r=-7.6998-0.5572+1.1335i-0.5572-1.1335i0.8141例求多项式x4+8x3+3x2+4x-10的根P=[1,8,3,4,-10];r=roots(P)多项式根的r向量：[r(1),r(2),…,r(n)]（2）根行向量创建多项式（poly）P=poly(r)>>PP=poly(r)PP=1.00008.00003.00004.0000-10.0000>>formatrat避免浮点显示>>p=poly(r)p=1834-10多项式的向量系数：(3)多项式求值（polyval）Y=polyval(P,x)若x数值，则求多项式在该点x的值；若X是向量或矩阵，每元素x(i,j)多项式求值。代数多项式P求值数组乘方运算:Y=a0*X.^n+a1*X.^(n-1)…an+1>>p=[1,3,0,2];>>poly2str(p,'x')ans=x^3+3x^2+2>>a=[1,2;3,4]a=1234>>y=polyval(p,a)y=62256114>>y=polyval(p,2)y=22例例

Y=polyvalm(P,X)矩阵多项式求值（polyvalm）Y=a0*X^n+a1*X^n-1…an+1X必为方阵，作自变量代入多项式求值；结果阶数还是保持方阵阶数。相当于矩阵X乘方运算：设A为方阵，P代表多项式x3-5x2-2：pp=polyval(P,A)的含义pp=A.*A.*A

-5*A.*A-2Pm=polyvalm(P,A)的含义：Pm=A*A*A-5*A*A-2>>p=[1,-5,0,-2];>>a=[2,4;1,0]a=2410>>pp=polyval(p,a)pp=-14-18-6-2>>pm=polyvalm(p,a)pm=-18-8-2-14矩阵乘方运算:数组乘方运算:例例14(4)多项式加减运算相同次数多项式，加减其系数向量,

不同次数,低次多项式中高次项用0补足。

p2=0x3-

0x2+2x

+1

p1+p2=2x3-

x2+2x

+4

[2,-1,0,3]

[2,1]

[0,0,

2,1]

[2,-1,2,4]

p1=2x3-

x2+3例例(5)多项式乘法-向量卷积w=conv(P1,P2)多项式x4+8x3-10与2x2-x+3的乘积。>>p1=[1,8,0,0,-10];>>p2=[2,-1,3];>>p12=conv(p1,p2)p12=Columns1through5215-524-20Columns6through710-30例例>>h=[3,2,1,-2,1,0,-4,0,3];>>x=[1,-2,3,-4,3,2,1];>>y=conv(h,x);>>n=0:14;>>stem(n,y);%杆图>>xlabel('Timeindexn');%标坐标轴>>ylabel('Amplitude');>>title('OutputObtainedbyConvolution')%图形标题grid;卷积是两个变量在某范围内相乘后求和的结果。卷积的变量是序列x(n)和h(n)，y=conv(x,h)来实现卷级的，输出的结果个数等于x的长度与h的长度之和减去1。卷积公式：z(n)=x(n)*y(n)=∫x(m)y(n-m)dm.例(6)多项式除法-向量解卷积[Q,r]=deconv(P1,P2)Q：多项式P1除以P2的商式，r：P1除以P2的余式，Q和r仍是多项式系数向量。P1=conv(P2,Q)+rdeconv是conv的逆函数求多项式x4+8x3-10除以2x2-x+3的结果。>>[q,r]=deconv(p12,p2)q=1800-10r=Columns1through500000Columns6through700>>[q,r]=deconv(p2,p12)q=0r=2-13例(7)多项式的一阶导数（polyder）

k=polyder(P)k=polyder(P,Q)多项式乘积P·Q的导函数：多项式P的导函数：导函数的分子存入p，分母存入q。[p,q]=polyder(P,Q)多项式除法P/Q的导函数：例：求有理分式的导数。命令如下：P=[1];Q=[1,0,5];[p,q]=polyder(P,Q)例例21(8)多项式积分(polyint)I=polyint([2,-1,0,3],5);例：

p(x)

=2x3-

x2+3求常数项5k=polyint(P,c)k=polyint(P)不定积分，常数项取c不定积分，常数项取零例

4.2曲线拟合和数据插值4.2.1曲线拟合polyfit

4.2.2数据插值interp

拟合成多项式推算未知点值4.2.1.多项式曲线拟合（polyfit）

最小二乘法(误差平方和最小)；x、y已知数据的横、纵坐标；n为多项式次数;p=polyfit(x,y,n)1.N次多项式曲线拟合S结构体数组(struct)，估计预测误差，含R，df和normr。

R：先输入x构建范德蒙矩阵V，后QR分解，得上三角矩阵。

df：自由度，

df=length(y)-(n+1)。

df>0时，超定方程组求解，拟合点数比未知数(p(1)~p(n+1))多。

normr：标准偏差、残差范数，normr=norm(y-V*p)，

此处的p为求解之后的数值。

[p,S]=polyfit(x,y,n)年份1900191019201930194019501960197019801990人口75.9991.97105.71123.20131.66150.69179.32203.21226.50249.63对不同年份人口数据分别进行1次、2次和9次多项式拟合（polyfit），用poly2str表示多项式完整形式法；1次、2次和9次多项式估计2000年人口(polyval)，结合美国2000年人口普查截至2000年4月1日美国人口2.81421906亿数据绘制1900~2000年间的时间-人口数(polyval)曲线，要求用plot不同线型（LineSpec）绘制原始数据点、拟合的1次、2次和9次多项式曲线，标注图例，说明是否阶数越到高越好美国从1900年~1990年每隔10年进行人口普查的数据如下表所示：（百万）例线型说明数据点标记符说明颜色说明-实线(默认)+加号符r红色--双划线*星号g绿色-.点划线.实心圆b蓝色:虚线o空心圆c青绿色x叉号符m洋红色

s正方形y黄色

d菱形k黑色

^上三角形w白色

v下三角形

>

右三角形

<

左三角形

p五角星

h六边形

plot(x,y,‘LineSpec’)线型属性例如'r-.*'、'-.r*'、'*-.r'等形式是等效的属性控制符顺序不受限制，可有可无。plot(x,y,'-gh')clearn=1900:10:1990;r=[7599,9197,10571,12320,13166,15069,17932,20321,22650,24963];nrf1=polyfit(n,r,1)%1次多项式拟合nrfs1=poly2str(nrf1,'n')%1次多项式完整字符串表达式nrf2=polyfit(n,r,2)%2次多项式拟合nrf9=polyfit(n,r,9)%9次多项式拟合nrf10=polyfit(n,r,10)%10次多项式拟合nrf1_2000=polyval(nrf1,2000)%一次拟合得到的2000年人口数nrf2_2000=polyval(nrf2,2000)%2次拟合得到的2000年人口数nrf9_2000=polyval(nrf9,2000)%9次拟合得到的2000年人口数n20=1900:4:2000;%1900年到2000年的线性等分数组nrfv1=polyval(nrf1,n20);%从1900年到2000年间1次拟合人口数nrfv2=polyval(nrf2,n20);%从1900年到2000年间2次拟合人口数nrfv9=polyval(nrf9,n20);%从1900年到2000年间9次拟合人口数

plot(n,r,'or',n20,nrfv1,'-<b',n20,nrfv2,'-*k')holdonplot(n20,nrfv9,'-dg')legend('原始数据’,‘1次曲线','2次多项式曲线','9次多项式曲线');xlabel('年');ylabel('人口数')2009版多项式次数要适当，过低误差大，过高波动明显2014版Goodnessoffit

适合度

SSE拟合误差

（误差平方和）

RMSErootmeansquareerror均方根误差

Rsquare

方程的确定系数，0~1之间，越接近1，表明方程的变量对y的解释能力越强。

2.曲线拟合工具箱cftool

R-square=1，说明拟合的结果很好CustomEquations：用户自定义的函数类型

Exponential：指数逼近，2种：a*exp(b*x)、a*exp(b*x)+c*exp(d*x)

Fourier：傅立叶逼近，7种，基础型是a0+a1*cos(x*w)+b1*sin(x*w)

Gaussian：高斯逼近，8种，基础型是a1*exp(-((x-b1)/c1)^2)

Interpolant：插值逼近，4种，linear、nearestneighbor、cubicspline、shape-PreservingPolynomial：多形式逼近，9种，linear~、quadratic~、cubic~、4-9thdegree~

工具箱提供的拟合类型Power：幂逼近，有2种类型，a*x^b、a*x^b+c

Rational：有理数逼近，分子、分母共有的类型是linear~、quadratic~、cubic~、4-5thdegree~；此外，分子还包括constant型

SmoothingSpline：平滑逼近（翻译的不大恰当，不好意思）

SumofSinFunctions：正弦曲线逼近，有8种类型，基础型是a1*sin(b1*x+c1)

Weibull：只有一种，a*b*x^(b-1)*exp(-a*x^b)拟合工具cftool

x12345678910y1311122832457080104x=1:10;y=[1,3,11,12,28,32,45,70,80,104];%1次多项式拟合p1=polyfit(x,y,1)Ps1=poly2str(p1,'x')xx=1:0.05:10;y1=polyval(p1,xx);%5次多项式拟合p2=polyfit(x,y,5)Ps2=poly2str(p2,'x')y2=polyval(p2,xx);%9次多项式拟合p3=polyfit(x,y,9)Ps3=poly2str(p3,'x')y3=polyval(p3,xx);>>p4=polyfit(x,y,11)y4=polyval(p4,xx);Warning:Polynomialisnotunique;degree>=numberofdatapoints.例holdon;plot(x,y,'*');plot(xx,y1,'r-');plot(xx,y2,'g-');plot(xx,y3,'b-')多项式次数要适当，过低误差大，过高波动明显已知：10个实验数据，分别采用2次、5次、9次和11次拟合，选出最佳拟合次数3.函数拟合>>clear>>x=[0:5];>>y=[0.2097,0.3523,0.4339,0.5236,0.7590,0.8998];>>ly=log(y);>>p=polyfit(x,ly,1);b=p(1);la=p(2);a=exp(la);Xx=linspace(0,5,30);Yy=a*exp(b*Xx);plot(x,y,'ro',Xx,Yy)例例插值构造的函数必须通过已知数据点；拟合则不要求,只要均方差最小。4.数据拟合和插值相同点都需根据已知数据构造函数；可使用得到函数计算未知点的函数值。不同点4.2.2数据插值构造函数近似表达式的方法。常用多项式作插值函数,称多项式插值。多项式插值基本思想：高次多项式或分段的低次多项式为被插值函数f(x)的近似解析表达式。1.插值法根据已知输入/输出数据集和当前输入估计输出值2.插值函数多项式插值法：拉格朗日插值法

牛顿插值法

埃尔米特插值法分段插值法样条插值法等

3.一维插值函数y=f(x)一维插值原理：调用格式：

yi=interp1(x,y,xi,'method')已知X,Y两个等长向量，采样点和样本值；Xi是向量或标量，欲插值点；Yi是一个与Xi等长的插值结果。采用差值法估计美国的人口数量2000年人口及1900~1996年每隔4的人口数(interp1)；将上题原始数据点、2阶拟合曲线和插值曲线绘制在同一图，标注坐标轴为‘时间’和‘人口’（xlabel、ylabel）、图形标题（title）为‘插值和拟合’和图例legend年份1900191019201930194019501960197019801990人口75.9991.97105.71123.20131.66150.69179.32203.21226.50249.63例clearn=1900:10:1990;r=[7599,9197,10571,12320,13166,15069,17932,20321,22650,24963];n4=1900:4:1996;%1900~1996年每隔4的线性等分数组nrf2=polyfit(n,r,2)%2次多项式拟合nrfv2=polyval(nrf2,n4)nr4=interp1(n,r,n4,'spline')%插值1900~1996年每隔4的人口数nr2000=interp1(n,r,2000)%

线性插值2000年人口数nr2000=interp1(n,r,2000,'spline')%

插值2000年人口数plot(n,r,'or',n4,nrfv2,'-*k',n4,nr4,'-dg')%绘图legend('原始数据','2阶多项式拟合曲线','插值曲线')%图例

例分段线性插值linear：插值点样本值落在两相邻数据点的连线上.

(缺省).最邻近插值法nearest：两点间插值点对应值就是离两点最近那点值。

3次多项式插值cubic：已知数据求3次多项式，用多项式求插值。

3次样条插值spline：每分段内构造一3次多项式，插值函数除满足插值条件和节点处光滑条件，再按照样条函数插值。插值方法method

yy=spline(x,Y,xx)clearx=0:1:10;y=sin(x);xx=0:0.2:10;yy=interp1(x,y,xx)subplot(1,4,1)plot(x,y,'-*',xx,yy,'dr');subplot(1,4,2);y2=interp1(x,y,xx,'nearest');plot(x,y,'-*',xx,y2,'dr');subplot(1,4,3);y3=interp1(x,y,xx,'cubic');plot(x,y,'-*',xx,y3,'dr')subplot(1,4,4);y4=interp1(x,y,xx,'spline');plot(x,y,'-*',xx,y4,'dr')X1取值范围不能超出X给定范围，否则会给出“NaN”错误。例38.025 22138.075 26738.125 31838.175 48738.225 81538.275 150238.325 316238.375 620738.425 841438.475 815638.525 597538.575 347338.625 160138.675 71038.725 38638.775 32538.825 24738.875 23238.925 20738.975 18239.025 157插值解法XRD：nano-sic/Alcomposiste例plot(xrd1(:,1),xrd1(:,2))>>xx=linspace(38,39,80);yy=interp1(xrd1(:,1),xrd1(:,2),xx,'spline');>>plot(xx,yy,'-*g')>>x=xy(:,1)x=38.025038.075038.125038.175038.225038.275038.325038.375038.425038.475038.525038.575038.625038.675038.725038.775038.825038.875038.925038.975039.0250>>y=xy(:,2)y=22126731848781515023162620784148156597534731601710386325247232207182157yys=spline(xrd1(:,1),xrd1(:,2),xx)4.二维插值函数z=f(x,y)二维插值的原理：向量的采样点X,Y，与采样点对应函数值Z；向量或标量的欲插值点Xi,Yi，插值结果Zi。指定插值方法method：

‘linear’‘nearest’‘cubic’‘spline’

Xi,Yi取值不超X,Y给定范围，否则“NaN”错误。

zi=interp2(x,y,z,xi,yi,method)函数interp2二维插值：例运行结果如下图所示。向量的采样点X,Y，与采样点对应函数值Z；向量或标量的欲插值点Xi,Yi，插值结果Zi。指定插值方法method：linear(缺省)nearestcubicv4griddata算法

Xi,Yi取值不超X,Y给定范围，否则“NaN”错误。

[xi,yi,zi]=griddata(x,y,z,xi,yi,method)函数griddata二维栅格插值：输入参量（XI,YI）通常是规则的格点>>clear>>p=rand(30,3)

X=p(:,1)Y=p(:,2)Z=p(:,3)[Xi,Yi]=meshgrid(linspace(min(X),max(X),100))[Xi,Yi,Zi]=griddata(X,Y,Z,Xi,Yi,'cubic')mesh(Xi,Yi,Zi)[Xi,Yi,Zi]=griddata(X,Y,Z,Xi,Yi,'v4')mesh(Xi,Yi,Zi)

随机生成30*3矩阵例4.3函数极值1.函数最小值和零点调用格式一样[x,fval,exitflag,output]=

fminbnd(fun,x1,x2,options)最小值的解或零点根x、该点的函数值fval、程序退出标志exitflag输出结果output、待求根函数fun、初始值x0、左右边界x1,x2fval,exitflag,output和options可省略[x,fval,exitflag,output]=fminsearch(fun,x0,options)[x,fval,exitflag,output]=fzero(fun,x0,options)迭代解法输出量x：最小值的x解或零点的根(自变量值)；fval=f(x)：最小值或零点时函数值f(x)；fun：待求最小值或根函数f(x)：函数名(少用)；字符串表达式，匿名对象、函数句柄、内联函数；

exitflag：程序退出标志,

1收敛到解x；output：选定输出结果；x0：搜索初始值；x1，x2：左右边界

output=

迭代次数iterations:24函数评价次数funcCount:48算法algorithm:‘bisection,interpolation’对分插值退出信息message:'Zerofoundintheinterval[-28,190.51]'option：配置优化参数,

optimset函数定义参数;最优化工具箱的(20多个)选项设定;optimset命令:将option显示出来;改变其中某个选项;display:选项函数调用时中间结果显示方式‘off’不显示;‘iter’每步都显示;‘final’只显示最终结果。optimset(‘display’,‘off’)options=optimset('param1',value1...)[x,fval,exitflag,output]=

fminbnd(fun,x1,x2,options)2.求一元函数最小值求区间[x1x2]内函数f(x)最小值时x：X：返回最小值的解;fval=F(x)，即最小值；fun：用于定义需求解函数;x1，x2：范围；option：最优化工具箱的选项设定x=2.5148f=-0.4993e=1o=iterations:13funcCount:14algorithm:'goldensectionsearch,parabolicinterpolation'message:'优化已终止:当前的x满足使用1.000000e-04的OPTIONS.Tol...'functionf=mmfun(x)f=exp(-0.1*x)*(sin(x)^2)-0.5*(x+0.1)*sin(x)end[x,f,e,o]=fminbnd('mmfun',-10,10,optimset('display','iter'))求在（-10，10）最小值例[x,f,e,o]=fminbnd('mmfun',6,10,optimset('display','iter'))x=8.0236f=-3.5680e=1o=iterations:9funcCount:10algorithm:'goldensectionsearch,parabolicinterpolation'message:'优化已终止:当前的x满足使用1.000000e-04的OPTIONS.Tol...'黄金分割搜索，抛物线插值x=-10:0.05:10; y=exp(-0.1*x).*(sin(x).^2)-0.5*(x+0.1).*sin(x);plot(x,y)plot(x,y)Fminbnd在指定区域找到真解最小最小求在（6，10）最小值x=-10:0.05:10; y=exp(-0.1*x).*(sin(x).^2)-0.5*(x+0.1).*sin(x);plot(x,y)plot(x,y)[x,fval]=fminbnd('exp(-0.1*x)*(sin(x)^2)-0.5*(x+0.1)*sin(x)',-10,10)x=2.514797840754235fval=-0.499312445280039最小[x,fval]=fminbnd('exp(-0.1*x)*(sin(x)^2)-0.5*(x+0.1)*sin(x)',6,10)x=8.0236；fval=-3.5680Fminbnd在指定区域找到真解例[x,fval]=fminbnd('-exp(-0.1*x)*(sin(x)^2)+0.5*(x+0.1)*sin(x)',-8,-2)x=-4.830203748934343fval=-3.947274022275747最大值求解：-f(x)在区间(a,b)上的最小值x=-10:0.05:10; y=-exp(-0.1*x).*(sin(x).^2)+0.5*(x+0.1).*sin(x);plot(x,y)>plot(x,y)最大例3.求多元函数的最小值在初始x0附近寻找局部最小值；使用options选项来指定优化器的参数；fval,exitflag,output和options可以省略。[x,fval,exitflag,output]=fminsearch(fun,x0,options)著名Rosenbrock's"Banana"测试函数,理论极小值x=1,y=1.求极小值点

x=11fval=0exitflag=1iterations:24funcCount:49algorithm:'Nelder-Meadsimplexdirectsearch'message:'优化已终止:

当前的x满足使用1.000000e-04的OPTIONS.TolX的终止条件，...例[x,fval,exitflag,output]=fminsearch('fxy',[1;1])functionf=fxy(x)f=100*(x(2)-x(1).^2).^2+(1-x(1)).^2end[x,fval,exitflag,output]=fminsearch('100*(x(2)-x(1).^2).^2+(1-x(1)).^2',[1;1])方程代数方程微分方程线性方程非线性方程常微分方程偏微分方程4.5常微分方程的求解2.5线性方程组求解4.4非线性方程求解4.4非线性方程数值求解4.4.1单变量非线性方程求解(fzero)4.4.2非线性方程组的求解(fsolve)迭代解法4.4.1单变量非线性方程求解(fzero)[x,fval,exitflag,output]=fzero(fun,x1,x2,options)[x,fval,exitflag,output]=fzero(fun,x0,

options)2.调用格式：函数是否穿越横轴决定零点数值解，求方程f(x)=0根

1.解方程原理：fval,exitflag,output和options可省略3.求根函数fun【f(x)】的调用求f(x)=x-10x+2=0在x0=0.5附近的根。函数名：fzero（

‘funname’，x0）例最优化的结构体output，最优化取下列域：

algorithm使用算法

funcCount函数个数评估

interval

iterations发现区间的迭代次数

iterations发现零值点迭代次数

message退出信息

[x,f,e,o]=fzero('fun',0.5,optimset('display','iter'))x=0.3758；f=0；e=1o=intervaliterations:8iterations:5funcCount:21algorithm:‘bisection,interpolation‘二分法’message:'在区间[0.34,0.613137]中发现零'建立函数文件fun.m。functionfv=fun(x)fv=x-10.^x+2;x=-4:0.05:1;f=x-10.^x+2;plot(x,f)[x,y]=ginput(1)ans=-2.0012[x,fval]=fzero('fun',-1)x=-1.989761447718557fval=0(1)建立函数文件fun.m。functionfv=fun(x)fv=x-10.^x+2;(2)调用fzero函数求根。

[x,fval,e,output]=fzero('fun',0.5)x=0.3758fval=0>>[x,y]=fzero('fun',-3,0.9)x=-1.9898y=0例多个根只求x0最近的根[x,f,e,o]=fzero('fun',8,optimset('display','iter'))

x=NaNfval=NaNexitflag=-3正在退出fzero:终止搜索包含符号变化的区间因为在搜索期间遇到NaN或Inf函数值。请检查函数或使用其他起始值重试。o=intervaliter:22iterations:0funcCount:44algorithm:'bisection,interpolation'message:'正在退出fzero:终止搜索包含符号变化的区间因为在搜索期间遇到NaN或Inf函数值。例无根为Nan,函数句柄（Functionhandle)>>[x,y]=fzero(@fun,0.9)x=0.3758y=2.2204e-016匿名函数：F=@(变量表)函数表达式>>[x,fval]=fzero(@(x)x-10^x+2,-1,6)x=0.3758fval=2.2204e-016class(f)function_handlehfun=@+函数名fun：待求最小值或根函数f(x)：函数名(少用)；字符串表达式，匿名对象、函数句柄、内联函数；

字符串表达式：fzero(‘expression’,x0)>>[x,y,e,o]=fzero('x-10.^x+2',-1,6,optimset('Display','iter'))x=-1.9898y=0e=1o=intervaliterations:12iterations:6funcCount:31algorithm:'bisection,interpolation'message:'Zerofoundintheinterval[0.28,-2.28]'例例inline本质和function一样，它直接内嵌在命令行里，不用单独定义function。函数表达式内联函数inlineinline(‘expression’)>>f=inline('x-10.^x+2');>>f(3)ans=-995f=inline('x-10.^x+2','x')Z=fzero(f,0.5)或z=fzero(inline('x-10.^x+2'),0.5)hd=@funfa=@(x)x-10^x+2fs='x-10^x+2‘fi=inline('x-10^x+2')NameSizeBytesClassfa1x116function_handlefi1x1832inlinefs1x816charhd1x116function_handleclass(f)ans=inline例4.多项式求的根（roots）>>roots([1,0,-2,-5])ans=2.0946-1.0473+1.1359i-1.0473-1.1359i>>[x,fval]=fzero('x^3-2*x-5',3)x=2.0946fval=-8.8818e-016求解多项式例4.4.2非线性方程组求解x向量；F(x)

函数向量。x：返回的向量解;fval=F(x)：函数值向量；Jacobian：解x处的Jacobian阵；fun：用于定义需求解函数;x0：初始估计值，列向量；option：最优化工具箱的选项设定方程组的标准形式：F(x)=0[x,fval,exitflag,output,jacobian]=fsolve(fun,x0,options)functionF=myfun(x)F=[2*x(1)-x(2)-exp(-x(1));-x(1)+2*x(2)-exp(-x(2))];end[x,fval]=fsolve(@myfun,[-5;-5])求方程组解例FunctionF=myfun(x)。%x为自变量所构成的数组。F=[表达式1；表达式2；…表达式m]在x1=-5,x2=-5解

(1)建立函数文件myfun.m

functiony=myfun(x)y(1)=x(1)-0.6*sin(x(1))-0.3*cos(x(2));y(2)=x(2)-0.6*cos(x(1))+0.3*sin(x(2));[x,fval=fsolve('myfun',[0.5,0.5]',optimset('Display','iter'))在(0.5,0.5)附近解。例(2)初值x0=0.5,y0=0.5下，调用fsolve求方程的根。>>x=fsolve('myfun',[0.5,0.5]',optimset(‘display','iter'))x=0.63540.3734将求得的解代回原方程，检验结果是否正确，

可见得到了较高精度的结果。q=myfun([0.6354,0.3734])q=1.0e-004*-0.2744-0.6564例4.5常微分方程微分方程：y’=f(t,y),t0≤t≤tn初始条件：y(t0)=y01.一阶常微分方程：方程的初值问题数值解：求解y(t)在节点t0,t1,t2,t3….tn,处近似值y0,y1,y2，y3….yn；采用等距节点tn=t0+nh，n=0,1,..n,h是步长微分方程描述：(ODE)t：时间列向量;y:对应t的数值解(列向量);odefun：显式方程y’=f(t,y),或含混合矩阵方程M(t,y),tspan：时间t范围，形式t0,tf,或[t0,tf],

y0：初始条件的值,options：用odeset设置可选参数龙格－库塔法

[t,y]=ode23(odefun,tspan,y0,options)[t,y]=ode45(odefun,tspan,y0,options)一阶常微分方程数值解

设有初值问题，y’=(y^2-t-2)/4/(t+1);0≤t≤1;y0(t=0)=2试求其数值解，并与精确解(y(t)=sqrt(t+1)+1)比较。(1)建立函数文件funt.m。functionyp=funt(t,y)yp=(y^2-t-2)/4/(t+1);(2)求解微分方程。t0=0;tf=1;y0=2;[t,y]=ode23('funt',[t0,tf],y0);%求数值解y1=sqrt(t+1)+1;%求精确解t'y'y1'y为数值解，y1为精确值，显然两者近似。t=linspace(0,1,11)例求解著名的Rossler微分方程组a=b=0.2,c=5.7,x(0)=y(0)=z(0)=0functionddx=rossler(t,x,flag,a,b,c)ddx=[-x(2)-x(3);x(1)+a*x(2);b+(x(1)-c)*x(3)];end>>a=0.2;b=0.2;c=5.7;>>x0=[0,0,0]';[t,y]=ode45('rossler',[0,100],x0,[],a,b,c)>>plot(t,y)>>figure>>plot3(y(:,1),y(:,2),y(:,3))例例2.二阶常微分方程：x’’=f(t,x,x’)

求解振荡器的经典Verderpol的微分方程令y(1)=x,y(2)=dx/dty=[y(1),y(2)],y’=[y(1)’,y(2)’]一阶微分方程组:初始条件：x(t=0)=1，y(1)(t=0)=1x’(t=0)=0,y(2)(t=0)=0Y0=[1;0]globalMUMU=1[t,y]=ode23('verderpol1',[0,40],[1;0]);plot(t,y)functionyy=verderpol1(t,y)globalMUyy=[y(2);

MU*(1-y(1).^2).*y(2)-y(1)];MU=32.高于2阶常微分方程的求解基本过程

（1）规律、定律、公式的描述形式：初始条件：微分方程：（2）高阶方程(组)转换一阶：一阶微分方程组：

初始条件：

列向量（3）根据(1)与(2)，编写导数M函数文件；（4）将M文件与初始条件传递给求解器Solver;（5）运行后得到ODE的、在指定时间区间解列向量y（其中包含y及不同阶的导数）。3.显式常微分方程组解法对比[t,y]=solver(odefun,tspan,y0)solversolver求解器solver奇异矩阵奇异矩阵4.求解器Solver与方程组的关系函数指令含

义函

数含

义求解器Solverode23普通2-3阶法解ODEodefile包含ODE的文件ode23s低阶法解刚性ODE选项odeset创建、更改Solver选项ode23t解适度刚性ODEodeget读取Solver的设置值ode23tb低阶法解刚性ODE输出odeplotODE的时间序列图ode45普通4-5阶法解ODEodephas2ODE的二维相平面图ode15s变阶法解刚性ODEodephas3ODE的三维相平面图ode113普通变阶法解ODEodeprint在命令窗口输出结果5.不同求解器Solver的特点求解器SolverODE类型特点说明ode45非刚性一步算法；4，5阶Runge-Kutta方程；累计截断误差达(△x)3大部分场合的首选算法ode23非刚性一步算法；2，3阶Runge-Kutta方程；累计截断误差达(△x)3使用于精度较低的情形ode113非刚性多步法；Adams算法；高低精度均可到10-3～10-6计算时间比ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法；Gear’s反向数值微分；精度中等若ode45失效时，可尝试使用ode23s刚性一步法；2阶Rosebrock算法；低精度当精度较低时，计算时间比ode15s短ode23tb刚性梯形算法；低精度当精度较低时，计算时间比ode15s短6.常微分方程组解法器参数4.6数值积分和微分4.6.1数值积分

4.6.2数值微分将区间[a,b]等分n个子区间[xi,xi+1]，i=1…n，x1=a，xn+1=b。求定积分就求和问题。数值积分方法：简单的梯形法trapz

辛普生(Simpson)法

牛顿－柯特斯(Newton-Cotes)法4.6.1数值积分1.数值定积分基本原理向量自变量X和应变量Y(1)梯形法数值积分trapz用trapz函数计算定积分。

X=1:0.01:2.5;Y=exp(-X);%生成函数关系数据向量trapz(X,Y)ans=0.285796824163932.数值积分的实现方法z=trapz(Y)z=28.5797Z=trapz(Y)Z=trapz(X,Y)

y向量和>>[q,n]=quad('exp(-x)',1,2.5)q=0.2858n=13例

fun：被积函数；n：被积函数调用次数；a和b：定积分下限和上限；tol：控制积分精度，缺省1e-6；trace[q,n]=quad(fun,a,b,tol,trace)q=quad(fun,a,b,tol,trace)

(2).变步长辛普生法（自适应Simpleson）非0展现积分过程；0不展现，缺省时trace=0；返回参数I即定积分值是否展现积分过程P为压力kpa，V为体积，m3；n为摩尔数kmol；R为理想气体常数8.314kpam3/kmolK;T为温度。气缸内1mol气体，温度为300k，气体在整个过程恒温，体积由V1=1m3膨胀到V2=5m3求解气缸活塞做功n=input('n=');T=input('T=');R=8.314;pp=@(v)n*R*T./v;W=quad(pp,1,5)例

(1)建立被积函数fesin.m。functionf=fesin(x)f=exp(-0.5*x).*sin(x+pi/6);(2)调用数值积分函数quad

[S,n]=quad('fesin',0,3*pi)S=0.9008n=77

求定积分quad('(-0.5*x).*sin(x+pi/6)',0,3*pi)[q,n]=quad(‘sqrt(1+cos(x).^2)’,0,pi/2,1.0e-6,1))例

s=quad('0.2+sin(x)',0,pi/2,1e-8);>>x=0:pi/10:pi/2;>>y=0.2+sin(x);>>s1=sum(y);>>ss=s1*pi/10;>>st=trapz(x,y)formatshort>>disp(['quad积分',blanks(4),'sum求积分',blanks(4),'trapz求积分'])disp([s,ss,st])holdonplot(x,y,'linewidth',2)>>stairs(x,y,'-r*')>>stem(x,y,'-gh')quad积分sum求积分trapz求积分

1.31421.35781.3059例参数含义和quad相似；tol缺省值10-6；调用步数明显小于quad，高效率求值；积分精度更高。[I,n]=quadl(fun,a,b,tol,trace)(3)牛顿－柯特斯法

求定积分(1)被积函数文件fx.m。functionf=fx(x)f=x.*sin(x)./(1+cos(x).*cos(x));(2)调用函数quad8求定积分。I=quad8('fx',0,pi)I=2.4674

例分别用quad和quadl求定积分的近似值，并在相同的积分精度下，比较函数的调用次数。调用函数quad求定积分：formatlong;fx=inline('exp(-x)');[I,n]=quad(fx,1,2.5,1e-10)I=0.28579444254766n=65调用函数quad8求定积分：formatlong;fx=inline('exp(-x)');[I,n]=quad8(fx,1,2.5,1e-10)I=0.28579444254754n=33例在[xmin,xmax]×[ymin,ymax]区域二重定积分；参数tol，trace的用法与函数quad相同；tol或method可以忽略。3.二重定积分的数值求解q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method)

(1)建立函数文件fxy.m：functionf=fxy(x,y)globalki;ki=ki+1;%ki用于统计被积函数的调用次数f=exp(-x.^2/2).*sin(x.^2+y);(2)调用dblquad函数globalki;ki=0;I=dblquad('fxy',-2,2,-1,1)kiI=1.57449318974494ki=1038

计算二重积分f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y');dblquad(f,-2,2,-1,1)ans=1.5745dblquad(f,-2,2,-1,1)ans=1.5745dblquad('exp(-x.^2/2).*sin(x.^2+y)',-2,2,-1,1)例4.三重积分数值求解triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method)例计算三重积分>>triplequad(@(x,y,z)y.*sin(x)+z.*cos(x),0,pi,0,1,-1,1)ans=1.999999994362637>>triplequad(inline('y.*sin(x)+z.*cos(x)'),0,pi,0,1,-1,1)ans=1.999999994362637>>triplequad('y.*sin(x)+z.*cos(x)',0,pi,0,1,-1,1)tol或method可以忽略例4.7.2数值微分向前差分⊿f(x)=f(X+h)-f(x)向后差分⊿f(x)=f(X)-f(x-h)中心差分⊿f(x)=f(X+h/2)-f(x-h/2)h充分小差分1.数值差分与差商xX+hX-hX-h/2X+h/2f(x+h

）f(x）f(x-h）差商导数2.数值微分函数diff向量X向前差分：Y=diff(X)向量X的n阶向前差分：Y=diff(X,n)矩阵A的n阶差分：dX(i)=X(i+1)-X(i)，i=1,2,…,n-1。diff(X,2)=diff(diff(X))Y=diff(A,n,dim)dim=1时(缺省状态)，按行计算；dim=2按列计算。

生成以向量V=[1,2,3,4,5,6]为基础的范得蒙矩阵，按列进行差分运算。V=vander(1:6)DV=diff(V)%计算V的一阶差分V=1111116842181279312566416416251252551dv=15731065195101753771036961910例3.梯度函数gradient中心差分FX=gradient(F)[FX,FY]=gradient(F)一元函数：二元函数：Fx：向量F的中心梯度H：F相邻两点间的间距。F是二维矩阵，FX相当于dF/dx=dF/HXFY相当于dF/dy=dF/HYHX，HY各方向相邻点距离V=vander(1:5);[FX,FY]=gradient(V)FX=00000-8-6-3-3/2-1-54-36-12-4-2-192-120-30-15/2-3-500-300-60-12-4FY=1573104013410120286102724981036961910例clearde=5*eps;d=pi/100;x=0:d:2*pi;y=sin(x);yde=(sin(x+de)-sin(x))/de;yd=(sin(x+d)-sin(x))/d;ydd=diff(sin(x))/d;yg=gradient(sin(x))/d;holdonplot(x,y,x,yde,'-b*',x,yd,'yh',x(1:end-1),ydd,'r>',x,yg,'+k')已知y=sin(x),采用diff和gradient计算该函数在[0,2π]区间中的近似导函数。eps是一个极小的数:2.2204e-16，为了避免0带来运算错误用eps代替，用机器零阈值eps替代理论0计算极限例f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2');x=-3:0.01:3;p=polyfit(x,f(x),5);%用5次多项式p拟合f(x)dp=polyder(p);%对拟合多项式p求导数dpdpx=polyval(dp,x);%求dp在假设点的函数值f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2');dx=diff(f([x,3.01]))/0.01;%直接对f(x)求数值导数g=inline('(3*x.^2+4*x-1)./sqrt(x.^3+2*x.^2-x+12)/2+1/6./(x+5).^(5/6)+5');%

f(x)导数解析式gx=g(x);%求函数f的导函数g在假设点的导数plot(x,dpx,x,dx,'.',x,gx,'-');%作图

用不同方法求f(x)的导数，并在同一坐标中做f'(x)的图像。例4.7概率统计4.7.1二项分布的概率计算

4.7.2正态分布的概率计算

4.7.3专用概率函数列表

4.7.4数字特征

4.7.5交互界面统计

随机现象统计规律的数学方法1.二项分布

N次随机试验的随机现象满足条件：1)重复N次相互独立随机试验；2)每次试验两结果成功概率p失败的概率1-p。4.7.1二项分布bino的概率计算离散随机变量：重复N次试验取得p成功的次数为k，k可以在N+1个数即[0,1,…,n]范围取值的离散随机变量bino