MATLAB 语言程序设计基础.docx

《MATLAB 语言程序设计基础.docx》由会员分享，可在线阅读，更多相关《MATLAB 语言程序设计基础.docx（16页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、%第六章 微分方程问题的解法% 微分方程的解析解方法% 常微分方程问题的数值解法% 微分方程问题算法概述% 四阶定步长 Runge-Kutta算法及 MATLAB 实现% 一阶微分方程组的数值解% 微分方程转换% 特殊微分方程的数值解% 边值问题的计算机求解% 偏微分方程的解% 6.1 微分方程的解析解方法% y=dsolve(f1, f2, , fm ,x)% syms t; u=exp(-5*t)*cos(2*t+1)+5;% uu=5*diff(u,t,2)+4*diff(u,t)+2*u% syms t y;% y=dsolve(D4y+10*D3y+35*D2y+50*Dy+24*y

2、=.% 87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10,y(0)=3,Dy(0)=2,D2y(0)=0,D3y(0)=0)% x,y=dsolve(D2x+2*Dx=x+2*y-exp(-t),.% Dy=4*x+3*y+4*exp(-t)% syms t x; % x=dsolve(Dx=x*(1-x2)+1)% Warning: Explicit solution could not be found; implicit solution returned.% In D:MATLAB6p5toolboxsymbolicdsolve.m

3、at line 292% x =% t-Int(1/(a-a3+1),a=.x)+C1=0% 故只有部分非线性微分方程有解析解。% 6.2 微分方程问题的数值解法% 6.2.1 微方程问题算法概述%Euler算法%function outx,outy=MyEuler(fun,x0,xt,y0,PointNum)% fun 表示f(x,y); x0,xt:自变量的初值和终值; % y0:函数在x0处的值,其可以为向量形式; % PointNum表示自变量在x0,xt上取的点数% if nargin5 |PointNum=0 %PointNum 默认值为100% PointNum=100;% en

4、d% if nargin4 % y0默认值为0 % y0=0;% end% h=(xt-x0)/PointNum; %计算步长h% x=x0+0:PointNum*h; %自变量数组% y(1,:) = y0(:); %将输入存为行向量，输入为列向量形式% for k = 1:PointNum% f=feval(fun,x(k),y(k,:); f=f(:); %计算f(x,y)在每个迭代点的值 % y(k + 1,:) =y(k,:) +h*f; %对于所取的点x迭代计算y值% end% outy=y; outx=x; % plot(x,y) %画出方程解的函数图 % end% x1,y1=

5、MyEuler(myfun,0,2*pi,1,16);% myfun=inclineDy=sin(x)+y% function Xout,Yout=MyEulerPro(fun,x0,xt,y0,PointNumber)% %MyEulerPro 用改进的欧拉法解微分方程% if nargin5 | PointNumber=0 % %PointNumer默认值为100% PointNumer=100;% end% if nargin0，收敛于x，exitflag=0，% 表示超过函数估计值或迭代的最大数字，exitflag f=inline(exp(-2*t).*cos(10*t)+exp(-

6、3*(t+2)*sin(2*t),t); % 目标函数% t0=1; t1,f1=fminsearch(f,t0); t1 f1% ans =% 0.9228 -0.1547% t0=0.1; t2,f2=fminsearch(f,t0); t2 f2% ans =% 0.2945 -0.5436%7.2.4 利用梯度求解最优化问题 %7.2.5非线性最小二乘 % 函数 lsqnonlin% 格式 x = lsqnonlin(fun,x0) % x = lsqnonlin(fun,x0,lb,ub) % x = lsqnonlin(fun,x0,lb,ub,options)% %x0为初始解向

7、量；fun为 ，i=1,2,m， % lb、ub定义x的下界和上界, options为指定优化参数，若x没有界，则lb= ，ub= 。%由于lsqnonlin中的fun为向量形式而不是平方和形式，因此，myfun函数应由f(x)建立% function F = myfun10(x)% k = 1:10;% F = 2 + 2*k-exp(k*x(1)-exp(k*x(2);% 然后调用优化程序： % x0 = 0.3 0.4;% x,resnorm = lsqnonlin(myfun10,x0)% Optimization terminated: norm of the current ste

8、p is less% than OPTIONS.TolX.% x =% 0.2578 0.2578% resnorm =% 124.3622% % 7.3有约束最优化问题的计算机求解 % 7.3.1 约束条件与可行解区域 % 对于一般的一元问题和二元问题，可用图解法直接得出问题的最优解。 % 例：用图解方法求解：% x1,x2=meshgrid(-3:.1:3); % 生成网格型矩阵% z=-x1.2-x2; % 计算出矩阵上各点的高度% i=find(x1.2+x2.29); z(i)=NaN; % % 找出 x12+x229 的坐标，并置函数值为 NaN% i=find(x1+x21);

9、z(i)=NaN; % 找出 x1+x21的坐标,置为 NaN% surf(x1,x2,z); shading interp;% max(z(:) % 7.3.2 线性规划问题的计算机求解 %x,fopt,flag,c=linprog(f,A,B,Aeq,Beq,xm,Xm,x0,opt,p1,p2,.)% 7.3.3 二次型规划的求解%x,fopt,flag,c=quadprog(H,f,A,B,Aeq,Beq,xm,Xm,x0,opt,p1,p2,.)% 7.3.4 一般非线性规划问题的求解%x,fopt,flag,c=fmincon(F,x0,A,B,Aeq,Beq,xm,Xm,CF,o

10、pt,p1,p2,.)% 7.3.5 约束线性最小二乘 % 函数 lsqlin % 格式 x = lsqlin(C,d,A,b) % x = lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0,options) % %求在约束条件下，方程Cx = d的最小二乘解x。% Aeq、beq满足等式约束，若没有不等式约束，则设 A= ，b= 。% lb、ub满足，若没有等式约束，则Aeq= ，beq= 。% x0为初始解向量，若x没有界，则lb= ，ub= 。% options为指定优化参数% x,resnorm,residual,exitflag,output,lambda = lsql

11、in() % % resnorm=norm(C*x-d)2，即2-范数。% residual=C*x-d，即残差。exitflag为终止迭代的条件% output表示输出优化信. lambda为解x的Lagrange乘子% 7.4整数规划问题的计算机求解% 7.4.1 整数线性规划问题的求解% x,how=ipslv_mex(A,B,f,intlist,xM,xm,ctype)% A、B线性等式和不等式约束,，约束式子由ctype变量控制，% intlist为整数约束标示，how0表示结果最优，2为无可行解，其余失败。% % 7.4.2一般非线性整数规划问题与求解% err,f,x=bbnb2

12、sita(fun,x0,intlist,xm,xM,A,B,Aeq,Beq,CFun)% err字符串为空，则返回最优解。% 该函数尚有不完全之处，解往往不是精确整数，% 可用下面语句将其化成整数。% if (length(err)=0% x(intlist=1)=round(x(intlist=1)% end% 7.4.3 0-1规划问题求解% x=binprog(f,A,B,Aeq,Beq)% % 7.5极小化极大（Minmax）问题 % 函数 fminimax% 格式 x = fminimax(fun,x0)% x = fminimax(fun,x0,A,b,Aeq,beq,lb,ub,

13、nonlcon,options)% x,fval,maxfval,exitflag,output,lambda = fminimax()% 参数说明：nonlcon的作用是通过接受的向量x来计算非线性不等约束 ；% 先建立非线性约束函数，并保存为mycon.m：% function C,Ceq = mycon(x)% C = % 计算x处的非线性不等约束的函数值。% Ceq = % 计算x处的非线性等式约束的函数值。% fval为最优点处的目标函数值；% maxfval为目标函数在x处的最大值；% exitflag为终止迭代的条件；% lambda是Lagrange乘子，它体现哪一个约束有效。

14、% output输出优化信息。 % 第8章 概率论与数理统计问题的求解% 概率分布与伪随机数生成% 统计量分析% 数理统计分析方法及计算机实现% 统计假设检验% 方差分析及计算机求解% % 8.1 概率分布与伪随机数生成% 8.1.1 概率密度函数与分布函数概述% % 通用函数计算概率密度函数值 % 函数 pdf% 格式 P=pdf(name,K,A)% P=pdf(name,K,A,B)% P=pdf(name,K,A,B,C)% 说明 返回在X=K处、参数为A、B、C的概率密度值，对于不同的分布，参数个数是不同；name为分布函数名。% 例如二项分布：设一次试验，事件Y发生的概率为p，%

15、那么，在n次独立重复试验中，事件Y恰好发生K次的概率P_K为：% P_K=PX=K=pdf(bino,K,n,p) % pdf(norm,0.6578,0,1) %正态分布N（0，1）的随机变量X在点0.6578的密度函数值。% pdf(chi2,2.18,8) %自由度为8的卡方分布，在点2.18处的密度函数值。% % 随机变量的累积概率值(分布函数值) % 通用函数cdf用来计算随机变量的概率之和（累积概率值）% 函数 cdf% 格式 cdf(name,K,A)% cdf(name,K,A,B)% cdf(name,K,A,B,C)% % 说明 返回以name为分布、随机变量XK的概率之和

16、的累积概率值，name为分布函数名.% % cdf(norm,0.4,0,1) %标准正态分布随机变量X落在区间(-，0.4)内的概率。 % cdf(chi2,6.91,16) %求自由度为16的卡方分布随机变量落在0，6.91内的概率。 % % 随机变量的逆累积分布函数 % MATLAB中的逆累积分布函数是已知，求x。% 命令 icdf 计算逆累积分布函数% 格式 icdf(name,F,A)% icdf(name,F,A,B)% icdf(name,F,A,B,C) % 说明 返回分布为name，参数为a1,a2,a3，累积概率值为F的临界值，这里name与前面相同。% % 如果F= cd

17、f(name,X,A,B,C),% 则 X = icdf(name,F,A,B,C) % 在标准正态分布表中，若已知F=0.6554,求X% 解： icdf(norm,0.6554,0,1)% 例：公共汽车门的高度是按成年男子与车门顶碰头的机会% 不超过1%设计的。设男子身高X（单位：cm）服从正态分布N（175，6）% ，求车门的最低高度。% 解：设h为车门高度，X为身高。% 求满足条件 FXh=0.99，即 FX=0.01故% h=icdf(norm,0.99, 175, 6)% 8.1.2 常见分布的概率密度函数与分布函数 % 8.1.2.1 Poisson分布%y=poisspdf(x

18、,lamda)%F=poissCdf(x,lamda)%y=poissinv(F,lamda)% 例：绘制laml=1,2,5,10 时Poisson分布的概率密度函数与概率分布函数曲线。% x=0:15; y1=; y2=; lam1=1,2,5,10;% for i=1:length(lam1)% y1=y1,poisspdf(x,lam1(i); % y2=y2,poisscdf(x,lam1(i);% end% plot(x,y1), figure; plot(x,y2)%8.1.2.2 正态分布%y=normpdf(x,u,sigma)%F=normCdf(x,sigma)%y=no

19、rminv(F,sigma)% x=-5:.02:5; y1=; y2=;% mu1=-1,0,0,0,1; sig1=1,0.1,1,10,1; sig1=sqrt(sig1);% for i=1:length(mu1) % y1=y1,normpdf(x,mu1(i),sig1(i); % y2=y2,normcdf(x,mu1(i),sig1(i); % end% plot(x,y1), figure; plot(x,y2)% 8.1.2.3 T分布%y=gampdf(x,a,lamda)%绘制(alfa,laml)为（1,1),(1,0.5),(2,1),(1,2),(3,1)时% x

20、=-0.5:.02:5;x=-eps:-0.02:-0.5,0:0.02:5; x=sort(x);替代% y1=; y2=; % a1=1,1,2,1,3; lam1=1,0.5,1,2,1;% for i=1:length(a1)% y1=y1,gampdf(x,a1(i),lam1(i); % y2=y2,gamcdf(x,a1(i),lam1(i);% end% 8.1.2.4 卡方分布% y=chi2pdf(x,k)% F=chi2cdf(x,k) % x=chi2inv(F,k) % x=-eps:-0.02:-0.5,0:0.02:2; x=sort(x);% k1=1,2,3,

21、4,5; y1=; y2=;% for i=1:length(k1)% y1=y1,chi2pdf(x,k1(i);% y2=y2,chi2cdf(x,k1(i);% end% plot(x,y1), figure; plot(x,y2)% 8.1.2.5 t分布% y=tpdf(x,k)% F=tcdf(x,k) % x=tinv(F,k) % x=-5:0.02:5; k1=1,2,5,10; y1=; y2=;% for i=1:length(k1)% y1=y1,tpdf(x,k1(i); % y2=y2,tcdf(x,k1(i); % end% plot(x,y1), figure;

22、 plot(x,y2)% 8.1.2.6 Rayleigh分布% y=ray2pdf(x,b)% F=ray2cdf(x,b) % x=ray2inv(F,b) % x=-eps:-0.02:-0.5,0:0.02:1; x=sort(x);% p1=1 2 3 3 4; q1=1 1 1 2 1; y1=; y2=;% for i=1:length(p1)% y1=y1,fpdf(x,p1(i),q1(i);% y2=y2,fcdf(x,p1(i),q1(i); % end% plot(x,y1), figure; plot(x,y2)% 8.1.2.7 F 分布% y=f2pdf(x,a,

23、b)% F=f2cdf(x,a,b) % x=f2inv(F,a,b) % x=-eps:-0.02:-0.5,0:0.02:1; x=sort(x);% p1=1 2 3 3 4; q1=1 1 1 2 1; y1=; y2=;% for i=1:length(p1)% y1=y1,fpdf(x,p1(i),q1(i);% y2=y2,fcdf(x,p1(i),q1(i); % end% plot(x,y1), figure; plot(x,y2)% 8.1.3 概率问题的求解% % % 8.1.4 随机数与伪随机数% A=gamrnd(a,lamda,n,m) %生成n*m的T分布伪随机数

24、矩阵% B=chi2rnd(k,n,m) %生成x2分布伪随机数% C=trnd(k,n,m) %生成T分布伪随机数% D=frnd(p,q,n,m) %生成F分布伪随机数% E=raylrnd(b,n,m) %生成Rayleigh分布伪随机数% % 8.2 统计量分析 % 8.2.1 随机变量的均值与方差%定义% 求该向量各个元素的%均值m=mean、方差s2=var(x)和标准差s=std(x)、中位数m=median(x)% 在分布类型标识后加后缀 stat：% u,sigma2=gamstat(a,lamda)% 返回变量为相关分布的均值和方差%8.2.2 随机变量的矩%求给定随机向量

25、x的r阶原点矩与中心距：%A=sum（x.r)/length(x),B=moment(x,r)%8.2.3 多变量随机数的协方差分析%协方差矩阵%C=cov(x)%8.2.4 多变量正态分布的联合概率 密度即分布函数%p=mvnpdf(X,u,协方差矩阵）% 8.3数理统计分析方法及计算机实现 % 8.3.1 参数估计与区间估计% normfit(x,Pci)% 由极大拟然法估计出该分布的均值、方差 及其置信区间。置信度越大，% 得出的置信区间越小，即得出的结果越接近于真值。 % 还有gamfit(), raylfit(), poissfit() ，unifit()（均匀分布） 等参数估计函数

26、%8.3.2 多元线性回归与区间估计%a,aci=regress(y,X,aerfa)% errorbar()用图形绘制参数估计的置信区间。%8.3.3 非线性函数的最小二乘参数 估计与区间估计%a,r,J=nlinfit(x,y,Fun,a0) %最小二乘拟合%r为参数下的残差构成的向量。J为各个Jacobi行向量构成的矩阵。% 8.4 统计假设检验% 8.4.1 正态分布的均值假设检验% H,s,uci=ztest(X,u,sigma,aerfa)% H为假设检验的结论，当H0时表示不拒绝H0假设，否则表示拒绝该假设。% s为接受假设的概率值，为其均值的置信区间。% H,s,uci=tte

27、st(X,u,aerfa)% 若未知正态分布的标准差时，可用此函数。% 8.4.2 正态分布假设检验% 由随机样本判定分布是否为正态分布，可用下面两个假设算法的函数。% H,s=jbtest(X,aerfa) %Jarque-Bera检验% s为接受假设的概率值，s越接近于0，则可以拒绝是正态分布的原假设.%H,s=lillietest(X,aerfa) %Lilliefors 检验%8.4.3 其它分布的Kolmogorov-Smirnov 检验% 此函数(Kolmogorov-Smirnov 算法)可对任意已知分布函数进行有效的假设检验。% H,s=kstest(X,cdffun,aerf

28、a)% 其中cdffun为两列的值，第一列为自变量，第二列为对应的分布函数的值。% % 8.5方差分析及计算机求解 % 8.5.1 单因子方差分析% p,tab,stats=anova1(X)% X为需要分析的数据，每一列对应于随机分配的一个组的测试数据，% 返回概率p，tab为方差分析表 。stats为统计结果量，为结构变量，包括每组均值等。 % 8.5.2 双因子方差分析% p,tab,stats=anova2(X)% % 8.5.3 多因子方差分析% 采用manoval()函数进行多因子方差分析% % 8.6 统计作图% 8.6.1 正整数的频率表% 命令 正整数的频率表% 函数 tab

29、ulate% 格式 table = tabulate(X) % %X为正整数构成的向量，返回3列：% 第1列中包含X的值,% 第2列为这些值的个数，% 第3列为这些值的频率。% 8.6.2 经验累积分布函数图形% 函数 cdfplot% 格式 cdfplot(X) % %作样本X（向量）的累积分布函数图形% h=cdfplot(X) %h表示曲线的环柄图% h,stats=cdfplot(X) %stats表示样本的一些特征% % 8.6.3 最小二乘拟合直线% 函数 lsline% 格式 lsline %最小二乘拟合直线% h = lsline %h为直线的句柄% % 8.6.4 绘制正态分

30、布概率图形% 函数 normplot% 格式 normplot(X) % %若X为向量，则显示正态分布概率图形，若X为矩阵，则显示每一列的正态分布概率图形。% h = normplot(X) % %返回绘图直线的句柄% 说明 样本数据在图中用“+”显示；如果数据来自正态分布，则图形显示为直线，而其它分布可能在图中产生弯曲。% % 8.6.5 绘制威布尔(Weibull)概率图形% 函数 weibplot% 格式 weibplot(X) % %若X为向量，则显示威布尔(Weibull)概率图形，% %若X为矩阵，则显示每一列的威布尔概率图形。% h = weibplot(X) %返回绘图直线的柄% 说明绘制威布尔(Weibull)概率图形的目的是用图解法估计来自% 威布尔分布的数据X，如果X是威布尔分布数据，其图形是直线的，否则图形中可能产生弯曲% % 8.6.6 样本数据的盒图% 函数 boxplot% 格式 boxplot(X) % %产生矩阵X的每一列的盒图和“须”图，“须”是从盒的尾部延伸出来，并表示盒外数据长度的线，如果“须”的外面没有数据，则在“须”的底部有一个点。% boxplot(X,notch) % %当notch=1时，产生一凹盒图，notch=0时产生一矩箱图。% boxplot(X,notch,sym) % %sym