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*
2、y=.% 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 nargin0，收敛于 x，exitflag=0，% 表示超过函数估计值或迭代的最大数字，exitf
4、lag f=inline(exp(-2*t).*cos(10*t)+exp(-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 = ls
5、qnonlin(fun,x0,lb,ub,options)% %x0 为初始解向量；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(myfun
6、10,x0)% Optimization terminated: norm of the current step 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(
7、i)=NaN; % % 找出 x12+x229 的坐标，并置函数值为 NaN% i=find(x1+x21); 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 一般非线性规
8、划问题的求解%x,fopt,flag,c=fmincon(F,x0,A,B,Aeq,Beq,xm,Xm,CF,opt,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= 。
9、% options 为指定优化参数% x,resnorm,residual,exitflag,output,lambda = lsqlin() % % 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 变量控制
10、，% intlist 为整数约束标示，how0 表示结果最优，2 为无可行解，其余失败。% % 7.4.2 一般非线性整数规划问题与求解% err,f,x=bbnb2sita(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 极小化极大（ M
11、inmax）问题 % 函数 fminimax% 格式 x = fminimax(fun,x0)% x = fminimax(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)% x,fval,maxfval,exitflag,output,lambda = fminimax()% 参数说明：nonlcon 的作用是通过接受的向量 x 来计算非线性不等约束 ；% 先建立非线性约束函数，并保存为 mycon.m：% function C,Ceq = mycon(x)% C = % 计算 x 处的非线性不等约束的函数值。% Ceq = % 计算 x 处的非线性等式约束
12、的函数值。% fval 为最优点处的目标函数值；% maxfval 为目标函数在 x 处的最大值；% exitflag 为终止迭代的条件；% lambda 是 Lagrange 乘子，它体现哪一个约束有效。% output 输出优化信息。% 第 8 章 概率论与数理统计问题的求解% 概率分布与伪随机数生成% 统计量分析% 数理统计分析方法及计算机实现% 统计假设检验% 方差分析及计算机求解% % 8.1 概率分布与伪随机数生成% 8.1.1 概率密度函数与分布函数概述% % 通用函数计算概率密度函数值 % 函数 pdf% 格式 P=pdf(name,K,A)% P=pdf(name,K,A,B
13、)% P=pdf(name,K,A,B,C)% 说明 返回在 X=K 处、参数为 A、B、C 的概率密度值，对于不同的分布，参数个数是不同；name 为分布函数名。% 例如二项分布：设一次试验，事件 Y 发生的概率为 p，% 那么，在 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 处的密度函数值。% % 随机变量的累积
14、概率值( 分布函数值) % 通用函数 cdf 用来计算随机变量的概率之和（累积概率值）% 函数 cdf% 格式 cdf(name,K,A)% cdf(name,K,A,B)% cdf(name,K,A,B,C)% % 说明 返回以 name 为分布、随机变量 XK 的概率之和的累积概率值，name 为分布函数名.% % cdf(norm,0.4,0,1) %标准正态分布随机变量 X 落在区间(- ，0.4)内的概率。 % cdf(chi2,6.91,16) %求自由度为 16 的卡方分布随机变量落在0，6.91内的概率。 % % 随机变量的逆累积分布函数 % MATLAB 中的逆累积分布函数是
15、已知，求 x。% 命令 icdf 计算逆累积分布函数% 格式 icdf(name,F,A)% icdf(name,F,A,B)% icdf(name,F,A,B,C) % 说明 返回分布为 name，参数为 a1,a2,a3，累积概率值为 F 的临界值，这里 name与前面相同。% % 如果 F= cdf(name,X,A,B,C),% 则 X = icdf(name,F,A,B,C) % 在标准正态分布表中，若已知 F=0.6554,求 X% 解： icdf(norm,0.6554,0,1)% 例：公共汽车门的高度是按成年男子与车门顶碰头的机会% 不超过 1%设计的。设男子身高 X（单位：c
16、m）服从正态分布 N（175，6）% ，求车门的最低高度。% 解：设 h 为车门高度， X 为身高。% 求满足条件 FXh=0.01 故% h=icdf(norm,0.99, 175, 6)% 8.1.2 常见分布的概率密度函数与分布函数% 8.1.2.1 Poisson 分布%y=poisspdf(x,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(
17、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=norminv(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);
18、 % 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=-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
19、(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,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,
20、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; 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
21、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,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
22、,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 分布伪随机数矩阵% 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
23、统计量分析% 8.2.1 随机变量的均值与方差%定义% 求该向量各个元素的%均值 m=mean、方差 s2=var(x)和标准差 s=std(x)、中位数 m=median(x)% 在分布类型标识后加后缀 stat：% u,sigma2=gamstat(a,lamda)% 返回变量为相关分布的均值和方差%8.2.2 随机变量的矩%求给定随机向量 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,协
24、方差矩阵）% 8.3 数理统计分析方法及计算机实现% 8.3.1 参数估计与区间估计% normfit(x,Pci)% 由极大拟然法估计出该分布的均值、方差 及其置信区间。置信度越大，% 得出的置信区间越小，即得出的结果越接近于真值。 % 还有 gamfit(), raylfit(), poissfit() ，unifit()（均匀分布） 等参数估计函数%8.3.2 多元线性回归与区间估计%a,aci=regress(y,X,aerfa)% errorbar()用图形绘制参数估计的置信区间。%8.3.3 非线性函数的最小二乘参数 估计与区间估计%a,r,J=nlinfit(x,y,Fun,a0
25、) %最小二乘拟合%r 为参数下的残差构成的向量。 J 为各个 Jacobi 行向量构成的矩阵。% 8.4 统计假设检验% 8.4.1 正态分布的均值假设检验% H,s,uci=ztest(X,u,sigma,aerfa)% H 为假设检验的结论，当 H0 时表示不拒绝 H0 假设，否则表示拒绝该假设。% s 为接受假设的概率值，为其均值的置信区间。% H,s,uci=ttest(X,u,aerfa)% 若未知正态分布的标准差时，可用此函数。% 8.4.2 正态分布假设检验% 由随机样本判定分布是否为正态分布，可用下面两个假设算法的函数。% H,s=jbtest(X,aerfa) %Jarqu
26、e-Bera 检验% s 为接受假设的概率值，s 越接近于 0，则可以拒绝是正态分布的原假设 .%H,s=lillietest(X,aerfa) %Lilliefors 检验%8.4.3 其它分布的 Kolmogorov-Smirnov 检验% 此函数(Kolmogorov-Smirnov 算法) 可对任意已知分布函数进行有效的假设检验。% H,s=kstest(X,cdffun,aerfa)% 其中 cdffun 为两列的值，第一列为自变量，第二列为对应的分布函数的值。% % 8.5 方差分析及计算机求解% 8.5.1 单因子方差分析% p,tab,stats=anova1(X)% X 为需
27、要分析的数据，每一列对应于随机分配的一个组的测试数据，% 返回概率 p， tab 为方差分析表 。stats 为统计结果量，为结构变量，包括每组均值等。 % 8.5.2 双因子方差分析% p,tab,stats=anova2(X)% % 8.5.3 多因子方差分析% 采用 manoval()函数进行多因子方差分析% % 8.6 统计作图% 8.6.1 正整数的频率表% 命令 正整数的频率表% 函数 tabulate% 格式 table = tabulate(X) % %X 为正整数构成的向量，返回 3 列：% 第 1 列中包含 X 的值,% 第 2 列为这些值的个数，% 第 3 列为这些值的频
28、率。% 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 绘制正态分布概率图形% 函数 normplot% 格式 normplot(X) % %若 X 为向量，则显示正态分布概率图形，若 X 为矩阵，则显示每一列的正态分布概率图形
29、。% h = normplot(X) % %返回绘图直线的句柄% 说明 样本数据在图中用“+ ”显示；如果数据来自正态分布，则图形显示为直线，而其它分布可能在图中产生弯曲。% % 8.6.5 绘制威布尔(Weibull)概率图形% 函数 weibplot% 格式 weibplot(X) % %若 X 为向量，则显示威布尔(Weibull)概率图形，% %若 X 为矩阵，则显示每一列的威布尔概率图形。% h = weibplot(X) %返回绘图直线的柄% 说明绘制威布尔(Weibull) 概率图形的目的是用图解法估计来自% 威布尔分布的数据 X，如果 X 是威布尔分布数据，其图形是直线的，否则
30、图形中可能产生弯曲% % 8.6.6 样本数据的盒图% 函数 boxplot% 格式 boxplot(X) % %产生矩阵 X 的每一列的盒图和“须”图， “须”是从盒的尾部延伸出来，并表示盒外数据长度的线，如果“须”的外面没有数据，则在“须”的底部有一个点。% boxplot(X,notch) % %当 notch=1 时，产生一凹盒图，notch=0 时产生一矩箱图。% boxplot(X,notch,sym) % %sym 表示图形符号，默认值为“+” 。% boxplot(X,notch,sym,vert) % %当 vert=0 时，生成水平盒图，vert=1 时，生成竖直盒图（默认
31、值 vert=1） 。% boxplot(X,notch,sym,vert,whis) % %whis 定义 “须”图的长度，默认值为 1.5，若 whis=0 则 boxplot 函数通过绘制 sym 符号图来显示盒外的所有数据值。% % 8.6.7 给当前图形加一条参考线 % 函数 refline% 格式 refline(slope,intercept) % % slope 表示直线斜率，intercept 表示截距% % 8.6.8 在当前图形中加入一条多项式曲线% 函数 refcurve% 格式 h = refcurve(p) % %在图中加入一条多项式曲线，h 为曲线的环柄，p 为多
32、项式系数向量，% %p=p1,p2, p3,pn，其中 p1 为最高幂项系数。% % 8.6.9 样本的概率图形% 函数 capaplot% 格式 p = capaplot(data,specs) % %data 为所给样本数据，specs 指定范围，p 表示在指定范围内的概率。% 说明 该函数返回来自于估计分布的随机变量落在指定范围内的概率% % 8.6.10 附加有正态密度曲线的直方图% 函数 histfit% 格式 histfit(data) %data 为向量，返回直方图和正态曲线。% histfit(data,nbins) % nbins 指定 bar 的个数，% 缺省时为 data
33、 中数据个数的平方根。 % % 8.6.11 在指定的界线之间画正态密度曲线% 函数 normspec% 格式 p = normspec(specs,mu,sigma) % %specs 指定界线，mu,sigma 为正态分布的参数 p 为样本落在上、下界之间的概率。%例子%画直方图% x=normrnd(0.5,1.4,30000,1);% m=mean(x);s=std(x);% xx=-5:0.3:5;% yy=hist(x,xx);bar(xx,yy/length(x)/0.3); hold on% x0=-5:0.1:5; y0=normpdf(x0,0.5,1.4); plot(x
34、0,y0);%由 normfit() 函数可以直接求出置信区间% x=10.4,10.2,12,11.3,10.7,10.6,10.9,10.8,10.2,12.1;% m1,s1,ma,sa=normfit(x,0.05); m1, ma% %采用 T 检验函数即可判定是否接受均值为 mean(x) 的检验% x=10.4,10.2,12,11.3,10.7,10.6,10.9,10.8,10.2,12.1;% mean(x)% H,p,ci=ttest(x,mean(x),0.05)% x=9.78,9.17,10.06,10.14,9.43,10.6,10.59,9.98,10.16,1
35、0.09,9.91,10.36;% H,p,c,d=jbtest(x,0.05); % H% %经确认满足正态分布，所以用 normfit() 函数即可以求出方差及方差的置信区间% m1,s1,ma,sa=normfit(x,0.05); s1,sa%最小二乘拟合问题% x=1.027,1.319,1.204,0.684,0.984,0.864,0.795,0.753,1.058,0.914,1.011,0.926;% y=-8.8797,-5.9644,-7.1057,-8.6905,-9.2509,-9.9224,-9.8899,-9.6364,-8.5883, -9.7277,-9.023,-9.6605;% f=inline(a(1)*exp(-a(2)*x).*cos(a(3)*x+pi/3)+a(4)*exp(-a(5)*x).* cos(a(6)*x+pi/4),a,x);% c,ci=nlinfit(x,y,f,1;2;3;4;5;6)% x1,ii=sort(x); y1=y(ii);% y2=f(c,x1);% plot(x1,y1,x1,y2)
