《2015软件编程综合实习报告》.docx

2015软件编程综合实习报告姓名: xxxxx 学号: xxxxx 专业: 数学与应用数学 年级: 201x-2 中国石油大学（华东）理学院计算数学系 2015 年 二 学期目录一、ode45、ode15s求解ode方程2二、ode23求解延迟微分方程6三、bvp4函数求解边值问题6四、pdepe函数求解偏微分方程8五、单因素方差分析（anova1）12六、多因素方差分析（anova2）14七、回归分析181、一元多项式回归：182、多元线性回归：22八、聚类分析25九、判别分析29十、一维抛物线微分方程的求解及图形表述30十一、抛物线微分方程的求解及图形表述331.使用分数步长法332.产生动态图像35十二、有限元法37一、ode45、ode15s求解ode方程1、ode45 ode是Matlab专门用于解微分方程的功能函数。该求解器有变步长（variable-step）和定步长（fixed-step）两种类型。不同类型有着不同的求解器，其中ode45求解器属于变步长的一种，采用Runge-Kutta算法；和他采用相同算法的变步长求解器还有ode23。ode45表示采用四阶，五阶Runge-Kutta单步算法,截断误差为(x)3。解决的是Nonstiff(非刚性)常微分方程。ode45是解决数值解问题的首选方法，若长时间没结果，应该就是刚性的，可换用ode23试试。下面，我们在matlab中对ode45 解决常微分方程求解进行具体的操作。例：求解首先，在matlab中创建一个包含公式的函数刚性再在主窗口里写如下代码，我们就可以得到函数的解得图像options = odeset(RelTol,1e-4,AbsTol,1e-4 1e-4 1e-5);T,Y = ode45(rigid,0 12,0 1 1,options); plot(T,Y(:,1),-,T,Y(:,2),-.,T,Y(:,3),.)2、ode15s如果我们在求解方程的过程中，发现运用ode45 求解运行的非常缓慢，那么很有可能这个方程式刚性的，这是我们将用ode15s求解次方程，将会事半功倍。下面，我们就ode15s求解方程解疑具体实例。例：刚性系统的一个例子是由张弛振荡的范德波方程提供。极限周期有部分地方解决方案组件慢慢改变，问题是相当硬，具有非常急剧变化的区域交替它不是生硬。同ode45一样，在函数中，T,Y = ode15s(vdp1000,0 3000,2 0); plot(T,Y(:,1),-o)由图我们可以看到，在某一小区间，其值变化非常的迅速。二、ode23求解延迟微分方程 ode23 是解决低阶的非刚性的常微分方程，对于误差限或解决中等刚性问题。例：求解器solver中的龙格-库塔法求解应用实例1。分别采用二、三、四、五阶龙格-库塔方法求解以下方程，编写程序文件ex1208a.m：%ex1208a.m 用ode23 得到微分方程解并计算出该算法运行时间fun =inline(-3*y2+2*x.2+3*x,x,y); %用inline构造函数f(x,y)x,y=ode23(fun,0,1,1); %可得到x,y输出向量值ode23(fun,0,1,1),hold on %可得到输出的函数图结果如图12-5a所示。三、bvp4函数求解边值问题 bvp4c是用来求解常微分方程的边值问题。 它的语法为： sol = bvp4c(odefun,bcfun,solinit)sol = bvp4c(odefun,bcfun,solinit,options)solinit = bvpinit(x, yinit, params)参数： 例：边值问题可以有多个解决方案，并初步推测的目的之一是要表明你想要的解决方案。二阶微分方程 具有满足边界条件恰好两个解在此之前解决这个问题bvp4c，你必须写微分方程为两个一阶微分方程的系统这里以及，该系统具有所要求的形式function dydx = twoode(x,y) dydx = y(2) -abs(y(1);function res = twobc(ya,yb) res = ya(1) yb(1) + 2;solinit = bvpinit(linspace(0,4,5),1 0);sol = bvp4c(twoode,twobc,solinit);x = linspace(0,4);y = deval(sol,x);plot(x,y(1,:);四、pdepe函数求解偏微分方程 pdepe是解决一维抛物线，椭圆偏微分方程初，边值问题的函数，它的语法为： sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)sol,tsol,sole,te,ie=pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options) 参数：mA parameter corresponding to the symmetry of the problem. m can be slab = 0, cylindrical = 1, or spherical = 2.pdefunA handle to a function that defines the components of the PDE.icfunA handle to a function that defines the initial conditions.bcfunA handle to a function that defines the boundary conditions.xmeshA vector x0, x1, ., xn specifying the points at which a numerical solution is requested for every value in tspan. The elements of xmesh must satisfy x0 x1 . = 3.tspanA vector t0, t1, ., tf specifying the points at which a solution is requested for every value in xmesh. The elements of tspan must satisfy t0 t1 . = 3.optionsSome options of the underlying ODE solver are available in pdepe: RelTol, AbsTol, NormControl, InitialStep, MaxStep, and Events. In most cases, default values for these options provide satisfactory solutions. See odeset for details. 例：这个例子说明了简单的配方，计算和绘制一个偏微分方程的解。 ， 方程满足的初始条件为： 边界条件为： 我们很方便的是使用子函数放置在一个单一的M文件需要由pdepe的功能:function pdex1m = 0;x = linspace(0,1,20);t = linspace(0,2,5);sol = pdepe(m,pdex1pde,pdex1ic,pdex1bc,x,t);% Extract the first solution component as u.u = sol(:,:,1);% A surface plot is often a good way to study a solution.surf(x,t,u) title(Numerical solution computed with 20 mesh points.)xlabel(Distance x)ylabel(Time t)% A solution profile can also be illuminating.figureplot(x,u(end,:)title(Solution at t = 2)xlabel(Distance x)ylabel(u(x,2)% -function c,f,s = pdex1pde(x,t,u,DuDx)c = pi2;f = DuDx;s = 0;% -function u0 = pdex1ic(x)u0 = sin(pi*x);% -function pl,ql,pr,qr = pdex1bc(xl,ul,xr,ur,t)pl = ul;ql = 0;pr = pi * exp(-t);qr = 1;便可得到如下图形：五、单因素方差分析（anova1） 单因素方差分析ANOVA的目的是找出是否从几组数据具有共同的均值。也就是说，为了确定组是否在所测量的特性实际上是不同的。 单因素ANOVA是线性模型的一个简单的特殊情况。所述单向ANOVA形式的模型是： 例：输入一个矩阵 hogg = 24 14 11 7 19 15 7 9 7 24 21 12 7 4 19 27 17 13 7 15 33 14 12 12 10 23 16 18 18 20 p,tbl,stats = anova1(hogg);得到图像：您可以使用F统计量做假设检验，以找出是否细菌数是相同的。的anoval返回此假设检验的p值。在这种情况下，在p值为约0.0001，一个非常小的值。这是一个强烈的迹象表明，从不同的出货量细菌数是不一样的。 F统计极端作为观 察F将偶然在一万次只出现一次，当计数真正平等。通过的anoval返回的P值取决于对随机扰动假设ij模型方程。用于p值是正确的，这些干扰必须是独立的，正态分布的，并且具有恒定的方差。你可以得到一些图形保证该装置是通过观察箱形图中通过的anoval显示的第二图形窗口不同。然而，请注意，凹口被用于中位数的比较，而不是手段的比较。此屏幕上的更多信息，请参见箱线图。该函数返回表示统计显著性的p值，以及两个图，一个图在前面已经说过，为方差分析表，另一个为直观的表示各个变量的均值和均值置信区间的箱线六、多因素方差分析（anova2） MATLAB中双因素试验的方差分析用函数anova2实现。其调用格式如下：P=anova2(X, reps)它执行平衡的双因素试验的方差分析，来比较X中两个或多个列或行的均值。不同列的数据代表某一因素的差异，不同行的数据代表另一因素的差异。如果每行、列对有多于一个的观察点，则变量reps指出每一单元观察点的数目，每一单元包含reps行。 双因素方差分析的目的是找出是否从几组数据具有共同的均值。单向ANOVA和双向方差分析的区别在于基团在两因素方差分析有两类限定特征，而不是 一个。假设一家汽车公司有两个工厂，每个工厂生产相同的三种型号的汽车。这是合理的，如果要求在汽车的汽油里程从工厂变化到工厂以及从模型来建模。有两种预测，厂家和型号，来解释城市的差异。有可能是在里程的整体差，由于在工厂的生产方法的差异。有可能是在不同的模型（不考虑厂）的城市由于设计规格不同的差异。这些效应被称为添加剂。最后，一个工厂可能使高里程汽车在一个模型中（或许是因为一个优越的生产线），但无法从其他工厂其它机型不同。这种效应被称为交互。它不可能检测到的相互作用，除非有重复观测工厂和汽车模型的某种组合。双因素ANOVA是线性模型的一个特例。双向ANOVA形式的模型是 mileage = 33.3000 34.5000 37.4000 33.4000 34.8000 36.8000 32.9000 33.8000 37.6000 32.6000 33.4000 36.6000 32.5000 33.7000 37.0000 33.0000 33.9000 36.7000cars = 3;p,tbl,stats = anova2(mileage,cars);p = 0.0000 0.0039 0.8411您可以使用F统计量做假设的测试，以找出是否里程是一样的跨模型，工厂模型厂对（调整添加剂后的效果）。 anova2返回这些测试的p值。为模型的效果的p值是零至四个小数位。这是一个强烈的迹象表明，从里程一种模式的不同而不同。 F统计极端作为观察F将偶然少于一次在10000次发生，如果油耗是真正的平等，从模型到模型。如果你使用了multcompare函数来执行多重比较测试，你会发现，每对三种模式是显著不同。为工厂效果的p值是0.0039，这也是非常显著。这表明一个工厂出执行其他在它产生的汽车的汽油里程。所观察到的p值指示F统计极端作为观察F将偶然约五分之四的1000倍发生如果气体里程是从工厂到工厂真正相等。似乎没有被工厂和模型之间的任何相互作用。在p值，0.8411，意味着所观察到的结果很可能（84选自100倍）鉴于没有互动。通过anova2返回的p值依赖的假设随机扰动ijk模型方程。用于p值是正确的这些干扰必须是独立的，正态分布的，并且具有恒定的方差。此外，anova2要求数据进行平衡，而在这种情况下，意味着必须有模型和工厂的每一种组合的相同数量的汽车。下一节讨论，支持非平衡数据与任意数量的预测的功能。七、回归分析回归分析（Regression Analysis）是研究一个变量Y与其他若干变量X之间相关关系的一种数学工具。它是在一组试验或观测数据的基础上，寻找被随机性掩盖了的变量之间的依存关系。粗略地讲，可以理解为用一种确定的函数关系去近似代替比较复杂的相关关系。这个函数称为回归函数，在实际问题中称为经验公式。回归分析所研究的主要问题就是如何利用变量X、Y的观察值（样本），对回归函数进行统计推断，包括对它进行估计及检验与它有关的假设等。回归分析包含的内容广泛。此处将讨论一元多项式回归、多元线性回归、非线性回归以及逐步回归。1、一元多项式回归： 如果从数据的散点图上发现y与x呈现较明显的二次（或高次）函数关系，则可以选用多项式回归。polyfit是多项式系数估计函数，其中用参数n设定多项式的最高次幂。polyval与polyconf分别用于计算预测输出及置信区间。1用函数polyfit估计模型参数，其具体调用格式如下：p=polyfit(x,y,n) p返回多项式系数的估计值，n设定多项式的最高次数。x，y为对应数据点值。p,S=polyfit(x,y,n) S包含x的范德蒙矩阵的QR分解的R对角线、自由度以及残差，用来作为polyval函数的输入，从而估计误差。2输出预估值与残差的计算用函数polyval实现，其具体调用格式如下：Y=polyval(p,X) Y为X在估计参数为p的线性模型作用下的输出。Y, DELTA = polyval(p, X, S) p，S为polyfit的输出，DELTA为误差估计。在线性回归模型中，YDELTA以50%的概率包含函数在X处的真值。3模型预测的置信区间用polyconf实现，具体调用格式如下：Y, DELTA = polyconf(P, X, S) 返回预测的95%置信区间为Y DELTA；Y, DELTA = polyconf(P, X, S, alpha, value1) 通过value1设定显著性水平的值。4交互式画图工具polytool，调用格式如下：polytool(x, y, n)polytool(x, y, n, alpha)用n次多项式拟合y、x的值，默认值为1，并作出显著性水平作为alpha的置信区间，默认值为0.05。【例13-27】 回归分析综合实例。将17至29岁的运动员每两岁一组分为7组，每组两人测量其旋转定向能力，以考察年龄对这种运动能力的影响。先得到的一组数据如表13-24所示。试建立两者之间的关系。解：数据的散点图明显地呈现两端低中间高的形状，所以应拟合一条二次曲线。选用二次模型。通过在MATLAB中编写程序ex1327.m来实现，具体代码如下所示。%ex1327.m 多项式拟合%输入数据x0=17:2:29;x0=x0,x0;y0=20.48,25.13,26.15,30.0,26.1,20.3,19.35.24.35,28.11,26.3,31.4,26.92,25.7,21.3;%多项式系数拟合，得系数pp,s=polyfit(x0,y0,2); p%y的拟合值及预测值y的置信半径deltay,delta=polyconf(p,x0,s);y%交互式画图polytool(x0,y0,2)所得系数p以及x0的预测值如下： ex1327p = -0.2003 8.9782 -72.2150y = Columns 1 through 10 22.5243 26.0582 27.9896 28.3186 27.0450 24.1689 19.6904 22.5243 26.0582 27.9896 Columns 11 through 14 28.3186 27.0450 24.1689 19.6904polytool所得的交互式图形如图2、多元线性回归： 在MATLAB统计工具箱中使用函数regress实现多元线性回归。具体调用格式为：b = regress(Y, X)b, bint, r, rint, stats = regress(Y, X, alpha)例：多元线性回归综合实例。已知某湖泊八年来湖水中COD浓度实测值（y）与影响因素，如湖区工业产值（x1）、总人口数（x2）、捕鱼量（x3）、降水量（x4）的资料，建立污染物y的水质分析模型。解：通过在MATLAB中编写程序ex1328.m来实现，具体代码如下所示。%ex1328.m 多元线性回归%输入数据x1=1.376, 1.375, 1.387, 1.401, 1.412, 1.428, 1.445, 1.477;x2=0.450, 0.475, 0.485, 0.500, 0.535, 0.545, 0.550, 0.575;x3=2.170 ,2.554, 2.676, 2.713, 2.823, 3.088, 3.122, 3.262;x4=0.8922, 1.1610 ,0.5346, 0.9589, 1.0239, 1.0499, 1.1065, 1.1387;y=5.19, 5.30, 5.60,5.82,6.00, 6.06,6.45, 6.95;%保存数据（以数据文件.mat形式保存，便于以后调用）save ex1328_data x1 x2 x3 x4 y%取出数据load ex1328_data %执行回归命令x=ones(8,1),x1,x2,x3,x4;b,bint,r,rint,stats=regress(y,x);b %输出参数bin %各参数估计的置信区间stats %几个特殊统计量执行后，输出如下结果： ex1328b = -13.9849 13.1920 2.4228 0.0754 -0.1897bint = -26.0019 -1.9679 1.4130 24.9711 -14.2808 19.1264 -1.4859 1.6366 -0.9638 0.5844stats = 0.9846 47.9654 0.0047 0.0123所以，回归模型方程是： 此外，由stats的值可知 =0.9846，F=47.9654，P=0.0123。八、聚类分析Matlab提供了两种方法进行聚类分析。一种是利用 clusterdata函数对样本数据进行一次聚类，其缺点为可供用户选择的面较窄，不能更改距离的计算方法；另一种是分步聚类：（1）找到数据集合中变量两两之间的相似性和非相似性，用pdist函数计算变量之间的距离；（2）用 linkage函数定义变量之间的连接；（3）用 cophenetic函数评价聚类信息；（4）用cluster函数创建聚类。1Matlab中相关函数介绍1.1 pdist函数调用格式：Y=pdist(X,metric)说明：用 metric指定的方法计算 X 数据矩阵中对象之间的距离。X：一个mn的矩阵，它是由m个对象组成的数据集，每个对象的大小为n。metric取值如下：euclidean：欧氏距离（默认）；seuclidean：标准化欧氏距离；mahalanobis：马氏距离；cityblock：布洛克距离；minkowski：明可夫斯基距离；cosine：correlation： hamming：jaccard： chebychev：Chebychev距离。1.2 squareform函数 调用格式：Z=squareform(Y,.) 说明： 强制将距离矩阵从上三角形式转化为方阵形式，或从方阵形式转化为上三角形式。1.3 linkage函数调用格式：Z=linkage(Y,method)说 明：用method参数指定的算法计算系统聚类树。 Y：pdist函数返回的距离向量； method：可取值如下： single：最短距离法（默认）； complete：最长距离法；average：未加权平均距离法； weighted： 加权平均法；centroid： 质心距离法； median：加权质心距离法；ward：内平方距离法（最小方差算法）返回：Z为一个包含聚类树信息的（m-1）3的矩阵。1.4 dendrogram函数调用格式：H，T，=dendrogram(Z,p，)说明：生成只有顶部p个节点的冰柱图（谱系图）。1.5 cophenet函数调用格式：c=cophenetic(Z,Y)说明：利用pdist函数生成的Y和linkage函数生成的Z计算cophenet相关系数。1.6 cluster 函数调用格式：T=cluster(Z,)说明：根据linkage函数的输出Z 创建分类。1.7 clusterdata函数调用格式：T=clusterdata(X,)说明：根据数据创建分类。T=clusterdata(X,cutoff)与下面的一组命令等价：Y=pdist(X,euclid);Z=linkage(Y,single);T=cluster(Z,cutoff);2. Matlab程序2.1 一次聚类法X=11978 12.5 93.5 31908;57500 67.6 238.0 15900;T=clusterdata(X,0.9)2.2 分步聚类Step1 寻找变量之间的相似性用pdist函数计算相似矩阵，有多种方法可以计算距离，进行计算之前最好先将数据用zscore函数进行标准化。X2=zscore(X); %标准化数据Y2=pdist(X2); %计算距离Step2 定义变量之间的连接Z2=linkage(Y2);Step3 评价聚类信息 C2=cophenet(Z2,Y2); /0.94698Step4 创建聚类，并作出谱系图 T=cluster(Z2,6); H=dendrogram(Z2);分类结果：加拿大，中国，美国，澳大利亚，日本，印尼，巴西，前苏联剩余的为一类。九、判别分析判别分析使用训练数据来估计的预测变量的判别函数的参数。判别函数确定各阶级之间的预测变量空间的界限。所得分类器鉴别基于所述预测数据的类（的响应的分类级别）之间。统计工具箱功能分类进行判别分析。例：load fisheririsSL = meas(51:end,1);SW = meas(51:end,2);group = species(51:end);h1 = gscatter(SL,SW,group,rb,v,off);set(h1,LineWidth,2)legend(Fisher versicolor,Fisher virginica,. Location,NW)十、一维抛物线微分方程的求解及图形表述例题1 -sint 0x0 A=1 ux0,x=ux1,x=0 t0 u(x,0)=cos x 0x1解：使用向前差分格式 ujn+1-ujnw=uj+1n-2ujn-uj-1nh2+sintn j=1J-1,n=1N得出的紧凑格式为1111+r*-221-211-212-2)u1nu2nuJ+1n+f=u1n+1u2n+1uJ+1n+1可得到的程序为xa = 0;xb = 1;T = 1; % the time interval is 0,1N = 12800; tau = T/N; % time step sizet = (0:tau:T); a = 1;r = 0.5;J = 80;h = (xb-xa)/J;x = (xa:h:xb); e1 = ones(J+1,1);Matrix_I = spdiags(e1,0,(J+1),(J+1);Matrix_A = spdiags(e1,-2*e1,e1,-1,0,1,(J+1),(J+1);Matrix_A(1,2) = 2;Matrix_A(end,end-1) = 2; X = cos(pi*x);C = zeros(J+1,N+1); % collect all the solutionsC(:,1) = X;for n=1:N f = sin(t(n)*e1; f = tau*f; X = (Matrix_I + r