[理学]第二题论文.doc

2010高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： 我们的参赛报名号为（如果赛区设置报名号的话）： 第七组 所属学校（请填写完整的全名）： 三峡大学 参赛队员 (打印并签名) ：1. 尹环 2. 牛威 3. 王蒙 指导教师或指导教师组负责人 (打印并签名)： 指导老师组 日期： 2010 年7月 21 日赛区评阅编号（由赛区组委会评阅前进行编号）：2010高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：自动化车床管理论文摘要本文建立的是自动化车床单刀具连续加工零件工序定期检查和刀具更换的随机性优化模型。由于刀具在制造和运行过程中各种随机因素的影响，它的寿命是随机变量，我们对题中所给数据进行处理和理论分析，并用卡方拟合检验法，确定了刀具寿命服从的正态分布。通过分析，单个零件的期望损失最小，则可使总的期望损失达到最小，进而使工序得到最好的效益，因此我们把求解生产单个零件的期望损失费用最小作为目标函数。对问题一：有了使单个零件的期望损失费用达到最小这个目标函数，我们根据需要找出未发生故障而采取预防性更换的损失费用和刀具发生故障而带来的损失费用表达式，再找出实际更换刀具的期望周期长度，最后求解出最小单个零件的损失费用的目标函数，利用Matlab对其进行编程求解到如下结果：检查间隔=19件，换到周期=334件，平均每个零件的损失费用=4.452元/件。对问题二：在问题一的基础上，根据问题二的条件，相应改变未发生故障而采取预防性更换的损失费用和刀具发生故障而带来的损失费用表达式，再结合题目所给数据建立类似模型 同样以单个零件的损失费用为目标函数，利用Matlab对其进行编程求解到如下结果：检查间隔=29件，换到周期=319件，平均每个零件的损失费用=9.128元/件。对问题三：通过问题二我们发现当刀具发生故障时，其合格品概率达40%，相当高，容易引起工人的误判，所以我们可增加每次检查零件时的件数来增加可靠性。在这里我们选择每次通过检查两个零件，来减少误判的概率以使工序得到最好的效益，利用Matlab对其进行编程求解到如下结果：检查间隔=31件，换到周期=341件，平均每个零件的损失费用=5.975元/件。 此三种模型在合理分析如何使工序得到最好的效益的基础上，模型改进通过进一步减少工人误判的概率，更是对该种方案的优化。关键词：卡方拟合正态分布检验 动态随机模型 预防性换刀周期 1.问题重述自动化车床管理用自动化车床连续加工某种零件，此工序会因为刀具损坏或其他原因导致工序故障，该故障由刀具损坏引起的可能性为95%，其他原因引起的概率为5%。工序出现故障是随机的，假定生产任一零件时出现故障的概率均相等。工作人员只能通过检查零件来确定工序是否出现故障。现积累有100次刀具故障记录，故障出现时该刀具完成的零件数如附表。现计划在刀具加工一定件数后定期更换新刀具。生产工序的费用参数如下：故障时产出的零件损失费用 f=200元/件；进行检查的费用 t=10元/次；发现故障进行调节使恢复正常的平均费用 d=3000元/次(包括刀具费)；未发现故障时更换一把新刀具的费用 k=1000元/次。问题一：假定工序故障时产出的零件均为不合格品，正常时产出的零件均为合格品，为该工序系统设计检查间隔（生产多少零件检查一次）和刀具更换策略，使得该工序得到最好的效益。问题二：如果该工序正常时产出的零件不全是合格品，有2%为不合格品；而工序故障时产出的零件有40%为合格品，60%为不合格品。工序正常而误认有故障停机产生的损失费用为1500元/次。对该工序设计效益最好的检查间隔和刀具更换策略。问题三：在2)的情况,改进检查方式获得更高的效益。100次刀具故障记录(完成的零件数)4593626245425095844337488155056124524349826407425657065936809266531644877346084281153593844527552513781474388824538862659775859755649697515628954771609402960885610292837473677358638699634555570844166061064841204476545643392802466875397905816217245315125774964684995446457645583787656667632177153108512.模型假设1假定生产任一零件出现故障的机会均等，且相互独立；2假定我们所设方案为预防性换刀方案；3. 工作人员一经检查发现不合格零件，即认为工序出现故障；4. 我们认为各刀具的质量一样；5出现非刀具故障产生次品时，对后面生产不产生影响，且无需进行故障恢复。6，当每次检查零件件数为二时，以第二次检查为级准来判断零件是否故障。3，符号说明：在一个换刀周期内完成的零件件数；：生产的零件件数；：在一个检查间隔内完成的零件件数（包含正品和次品）；：在一个换刀周期内生产的合格零件件数；: 每次检查的费用t=10元/次；:故障时产生的零件损失费用；：发现故障进行调节是恢复正常的平均费用；：未发故障时更换一把新刀具的费用元/次；：未发生故障而采取预防性更换的损失费用（元）；：刀具发生故障而带来的损失费用（元）；：一个周期内的平均损失（元）；单个零件的损失费用（元）；：工序正常而误认有故障停机生产的损失费用a=1500元；：实际换刀周期的平均长度；：工序正常时产出的合格零件个数所占的百分比=；：工序故障时产出的合格零件个数所占的百分比=；:可靠度，刀具寿命大于生产件零件的概率：；：每次检查时检查的零件个数的件数：=2（件）。4问题分析首先我们通过对100次刀具故障记录数据的分析以及卡方拟合检验得知，刀具的寿命近似服从的正态分布，并以此为基础建立各问题的数学随机动态模型。4.1 问题一分析 由于刀具损坏和其他故障使工序出现故障，工序出现故障是随机的，工作人员通过检查零件来确定工序是否出现故障，并计划在刀具加工一定零件数后定期更换刀具。对于每一把刀具其可能加工的零件数都是相互独立的，呈现出一个随机的分布。题目要求我们设计效益最好的检查间隔和刀具更换策略。效益最好即为单个零件的期望损失费最少。总的损失费用可能由刀具在未发生故障而采取预防性损失费用（包括检查费用）、刀具发生故障而带来的损失费用（包括检查费用）这二部分组成。在明确了刀具寿命服从正态分布以后，我们以该分布为基础，模拟出刀具在未发生故障而采取预防性更换的期望损失费用，刀具发生故障而带来的期望损失费用和总的费用的函数表达式。其中，刀具在未发生故障而采取预防性更换的损失费用为：检查的次数乘以单次检查费用加上更换一把刀的费用；刀具发生故障而带来的损失费用为：检查的次数乘以单次检查费用加上故障进行调节使其恢复正常的平均费用再加上产生的故障零件所损失的费用；总的费用为：这两部分的费用分别乘以与之相对应的概率之和。接下来就是求出每个零件的损失费用的函数表达式。我们利用所得到的刀具寿命概率分布函数求得实际换刀周期的平均长度，这样，实际的单个零件的损失费用为：总的损失费用除以实际换刀周期的平均长度。我们以这个函数为目标函数进行Matlab分析，得到最优解。4.2 问题二分析 问题二的处理类似问题一的处理，不同的是题目假定工序正常时也有2%的次品生成，而在故障时也会有40%的合格品，因此组成各损失费用的两部分的求解就出现了不同。以此为基础列出刀具在未发生故障而采取预防性更换的损失费用，刀具发生故障而带来的损失费用和总的费用的函数表达式。其中刀具在未发生故障而采取预防性更换的损失费用为：零件合格所产生的费用（同问题一）乘以合格品所占的百分比再加上零件不合格所产生的费用（同问题一）乘以不合格品所占的百分比。刀具发生故障而带来的损失费用为：（同上）。总的费用为：这两部分的费用分别乘以与之相对应的概率之和。同样利用所得到的刀具寿命概率分布函数求得周期的平均长度，这样，单个零件的期望损失费用为：总的损失费用除以实际换刀周期的期望平均长度。我们以这个函数为目标函数进行Matlab分析，进而得到最优解。4.3 问题三分析问题三的处理类似问题一和问题二的处理，不同的是，我们将以前的每次检查的零件个数有所增加，因为，我们通过问题二发现，当工序故障时产出的零件有40%的合格品，60%的不合格品，这两种情况的概率都很高。容易引起工人的误判从而会增加工厂的费用。解决的办法是我们通过增加每次检查零件的件数，来进一步提高工人通过检查零件是否合格来判断工序是否故障的精确度。同样我们列出刀具在未发生故障而采取预防性更换的损失费用，刀具发生故障而带来的损失费用和总的费用的函数表达式。目标函数思路同问题二，以单个零件的期望损失费用为目标函数进行Matlab分析，进而得到最优解。5模型的建立5.1模型准备5.1.1正态性的验证首先我们根据题目中表格中的数据用Excel进行数据分析，发现这100次刀具故障时完成的零件数近似服从的正态分布。我们就大胆的假设这100次刀具故障数据近似服从正态分布，再用卡方拟合检验法来进一步验证。我们用来作为检验统计量（下面的证明原理详见概率论与数理统计第三版第241页243页，这里只给出验证计算）。我们假设刀具使用寿命近似服从的正态分布，：表示总体的分布函数；：表示总体Y的分布函数不是；若，则成立，即假设成立；若，则成立，即假设不成立；分析题中所给数据可以知道：Y的最小值为84，最大值为1153，我们可以把这110个数据从83.5到1153.5分成10组，每组间隔，画出频数直方图如下：频数直方图数出落在每个小区间的数据频数，算出，如下表所示：组限频数频率83.5-190.530.030.01455.758.52190.5-297.540.040.0430297.5-404.570.070.09699.695.06404.5-511.5160.160.164116.4115.6511.5-618.5250.250.213121.3129.33618.5-725.5200.200.199819.9820.02725.5-832.5130.130.143314.3311.79832.5-939.570.070.07837.836.26939.5-1046.530.030.03084.75.321046.5-1153.520.020.016210011100101. 9卡方分布检验正态分布表格其中计算得：1. 9通常我们取，则9.236故有，因此在水平0.1下我们接受，即总体的概率密度函数服从正态分布。由此我们得出刀具寿命服从正态分布。5.1.2参数估计 估计出该刀具寿命的平均值为600，方差为196.6，置信区间为561.0，639.0方差的0.95置信区间为172.6,228.4.综上检验我们接受原假设，由于刀具损失故障占95%起决定性作用。我们认为，整个工序无故障生产零件个数的分布近似为刀具无故障生产零件个数的分布，即5.2模型建立5.2.1问题一模型的建立 由问题一分析得到，我们采取以正态分布为基础来模拟整个生产过程。我们需求出未发生故障而采取预防性更换的损失费用，刀具发生故障而带来的损失费用，周期的平均长度，最终求得单个零件的平均损失费用。其中：刀具在未发生故障而采取预防性更换的损失费用为：检查的次数乘以单次检查费用加上更换一把刀的费用；刀具发生故障而带来的损失费用为：检查的次数乘以单次检查费用加上故障进行调节使其恢复正常的平均费用再加上产生的故障零件所损失的费用；一个周期内的平均损失费用为：这两部分的费用分别乘以与之相对应的概率之和。利用所得到的刀具寿命概率分布函数求得周期的平均长度，这样，单个零件的损失费用为：总的损失费用期望除以周期的平均长度。图一刀具未发生故障而采取预防性更换的损失费用为：刀具发生故障的损失费用为：刀具更换周期的平均长度为：一个周期内的平均损失费用为：单个零件的平均损失费用为：我们的目标函数就是使达到最小，即目标函数为 5.2.2问题二模型的建立由问题二分析得到，我们需求出未发生故障而采取预防性更换的损失费用，刀具发生故障而带来的损失费用，周期的平均长度，最终求得单个零件的平均损失费用。其中：刀具在未发生故障而采取预防性更换的损失费用为：零件合格所产生的费用（同问题一）乘以合格品所占的百分比再加上零件不合格所产生的费用（同问题一）乘以不合格品所占的百分比。刀具发生故障而带来的损失费用为：（同上）。总的费用为：这两部分的费用分别乘以与之相对应的概率之和。同样利用所得到的刀具寿命概率分布函数求得周期的平均长度，这样，每个零件的损失费用为：总的损失费用除以周期的平均长度。刀具未发生故障而采取预防性更换的损失费用：刀具发生故障的损失费用：刀具更换周期的平均长度：一个周期内的平均损失：单个零件的平均损失费用：目标函数为 5.2.3问题三模型的建立在问题二的建模基础上，模拟思路类似，只是将检查次数改为二，每次的概率有所不同，具体建模过程为总的费用为：这两部分的费用分别乘以与之相对应的概率之和。同样利用所得到的刀具寿命概率分布函数求得周期的平均长度，这样，每个零件的损失费用为：总的损失费用除以周期的平均长度。刀具未发生故障而采取预防性更换的损失费用为：刀具发生故障的损失费用： 刀具更换周期的平均长度：一个周期内的平均损失：单个零件的平均损失费用：目标函数为 6 问题的解答及结果分析6.1问题一的解答及分析对问题一建立的模型，进行简单的求导后推出 取极小值的应满足： 我们用Matlab进行积分运算，得到以下结果：一个检查间隔内完成的零件件数，单位（件）在一个换刀周期内完成的零件件数，单位（件）单个零件的损失费用单位（元）193424.454 从表中数据我们得知当，当设计方案为：检查间隔=19件，即生产19件零件检查一次；换到周期=334件，即每生产334件零件换刀；单个零件的损失费用=4.454元。6.2问题二的解答及分析对问题二建立的模型，我们用Matlab进行积分运算，得到以下结果：一个检查间隔内完成的零件件数，单位（件）在一个换刀周期内完成的零件件数，单位（件）单个零件的损失费用单位（元）293199.128 从表中数据我们得知当，当设计方案为：检查间隔=29件，即生产29件零件检查一次；换到周期=319件，即每生产319件零件换刀；单个零件的损失费用=9.128元。问题一和问题二进行比较我们发现，6.3问题三的解答及分析对问题三建立的模型，我们用Matlab进行积分运算，得到以下结果：一个检查间隔内完成的零件件数，单位（件）在一个换刀周期内完成的零件件数，单位（件）单个零件的损失费用单位（元）313415.975 从表中数据我们得知当，当设计方案为：检查间隔=31件，即生产31件零件检查一次；换到周期=341件，即每生产341件零件换刀；单个零件的损失费用=5.975元。比较问题二结果和问题三，我们发现，改进方案后单个零件损失费用有明显的减小，有=9.128元减小为=5.975元。这样采取模型三工序效益进一步提高，因此模型三有明显的优势。7 模型的评价、改进及推广7.1模型评价 优点：（1） 本文所建模型思路明确，易于理解，具有可操作性，具有广泛的应用价值。（2） 我们的所建模型对题目中数据进行了卡方拟合正态性检验，使我们所建模型利用到的数据准确，可靠，最优。（3） 本文用到的数学方法一般都是概率统计和对函数求导，通过计算机模拟所得，比较简单，具有普遍性。（4） 本文所建模型也运用于不同的变量中，即变量改变照样能模拟得到对应的最优解，对以后工序长期生产有指导性的价值。 缺点： （1）本模型主要考虑了刀具故障，其他故障未作深入研究；（2）有些假设可能存在不合理性。7.2 模型的改进 由于第三问的模型假设为当每次检查零件件数为二时，以第二次检查为级准来判断零件是否故障，这样尽管误判的概率会减小，但我们仍是每次检查两个零件，误判的概率还是稍微大一些，所以我们可以增加每次检查的零件数来进一步增加检查的精确性，可以是每次检查三个或四个甚至更高来减少误判的概率。这样对工厂的效益会进一步增强。从而使该种方案得到优化。7.3模型推广 由于我们所建模型为动态随机模型，可以运用于不同的变量中，所以此模型可以运用于多个行业领域，例如各种机械零件更换，汽车轮胎更换领域等。8.参考文献1数学模型第二版 姜启源 高等教育出版社 19962概率论与数理统计第三版 盛骤 谢式千 潘承毅 高等教育出版社 2003 3数学建模与实验宋来忠 王志明 北京科学出版社 2005附录问题一程序：syms x u m c; c1=u./x.*10+1000; c2=u./x.*10+3000+(u-m).*200; p1=0.00203*exp(-0.0000129.*(x-600).2); y1=int(p1,x,0,u); c3=c2.*y1+c1.*(1-y1); y2=int(u.*p1,x,u,inf); y3=int(x.*p1,x,0,u); y=y2+y3; c4=c3./y c4 =(10*u/x+3000+200*u-200*m)*(1879048192/820562063968459375*pi(1/2)*15229631907254606(1/2)*erf(1/34359738368*15229631907254606(1/2)*u-75/4294967296*15229631907254606(1/2)+1879048192/820562063968459375*erf(75/4294967296*15229631907254606(1/2)*pi(1/2)*15229631907254606(1/2)+(10*u/x+1000)*(1-1879048192/820562063968459375*pi(1/2)*15229631907254606(1/2)*erf(1/34359738368*15229631907254606(1/2)*u-75/4294967296*15229631907254606(1/2)-1879048192/820562063968459375*erf(75/4294967296*15229631907254606(1/2)*pi(1/2)*15229631907254606(1/2)/(1879048192/820562063968459375*u*pi(1/2)*15229631907254606(1/2)-1879048192/820562063968459375*erf(1/34359738368*15229631907254606(1/2)*u-75/4294967296*15229631907254606(1/2)*u*pi(1/2)*15229631907254606(1/2)-64563604257983430656/820562063968459375*exp(-42833339739153579375/9223372036854775808-7614815953627303/590295810358705651712*u2+571111196522047725/36893488147419103232*u)+45097156608/32822482558738375*pi(1/2)*15229631907254606(1/2)*erf(1/34359738368*15229631907254606(1/2)*u-75/4294967296*15229631907254606(1/2)+64563604257983430656/820562063968459375*exp(-42833339739153579375/9223372036854775808)+45097156608/32822482558738375*erf(75/4294967296*15229631907254606(1/2)*pi(1/2)*15229631907254606(1/2) vpa(c4) ans = (10.*u/x+3000.+200.*u-200.*m)*(.50089433903710877308927024420163*erf(.35916569992135941929847473336317e-2*u-2.1549941995281565157908484001790)+.49973897687681912453571729041131)+(10.*u/x+1000.)*(.50026102312318087546428270958869-.50089433903710877308927024420163*erf(.35916569992135941929847473336317e-2*u-2.1549941995281565157908484001790)/(.50089433903710877308927024420163*u-.50089433903710877308927024420163*erf(.35916569992135941929847473336317e-2*u-2.1549941995281565157908484001790)*u-78.682170542635657954656940854112*exp(-4.6440000000000000566656097045204-.12900000000000000157404471401446e-4*u2+.15480000000000000188885365681735e-1*u)+300.53660342226526385356214652098*erf(.35916569992135941929847473336317e-2*u-2.1549941995281565157908484001790)+300.60024125443500777024257422268)问题一结果（部分）：4.485201982922814.482287693074734.482280616080424.482273550890954.482266497476794.4795016247974.479494546785494.479487480580514.479480426152534.476864424514814.476857345478294.476850278254.476843222800434.47437555729014.474368477220714.474361408961274.474354352482274.472034499174864.47202741806474.472020348766224.472013291249924.469840736993744.469833654834864.469826584489424.469819525927894.467793768130594.4677866849154.467779613514614.467772553899894.465893100318684.465886016038344.465878943574984.465871882899074.464138251434764.464131166081594.464124092547194.464117030802024.462528749296624.462521662862494.462514588248934.46250752542644.461064131464214.461057043940944.461049968240074.461042904332024.459743945044264.459736856423644.459729779627234.459722714625484.458567746498174.458560656771924.458553578871734.458546512768034.457535101453264.457528010613064.457520931600784.457513864386864.456645584517184.456638492554674.456631412421964.456624344089474.455898779095524.455891686002294.455884604740754.455877535281324.455294277212414.455287182980024.455280100581214.45527302998644.454831679334194.454824583954134.454817500409564.45481042867094.454510594195934.454503497659644.454496412960784.454489340069754.454330638630854.454323540929744.4543164550684.454309381016034.454291437402564.454284338527994.454277251494754.454270176273234.454392623040034.45438552298334.454378434769874.454371358370134.454633835675124.45462673442754.454619645025164.45461256743854.455014722882914.455007620435614.45500052983564.454993451053254.455534939524454.455527835868634.455520744062124.455513664075284.456194147592034.456187042718824.456179949696944.456172868496764.456992016056964.456984909957434.456977815711264.456970733288854.457928220719664.457921113384814.457914017905414.45790693425184.459002444062114.458995335482934.458988238761274.458981153867484.460214375102644.460207265270064.460200167297094.460193081154074.46156370925294.46155659815784.46154949892444.461542411523074.463050148177054.463043035810244.463035935307274.463028846638484.464673399653074.464666286005354.464659184223614.464652094278164.466433177436224.466426062498324.466418959428544.466411868197214.468329201124484.468322084887084.46831498051997问题二程序：syms x u m c; c1=0.98.*(u./x.*10+1000)+0.02.*(u./x.*10+200.*x+3000+1500); c2=0.4.*(u./x.*10+1000)+0.6.*(u./x.*10+3000+(u-m).*200); p1=0.00203*exp(-0.0000129.*(x-600).2); y1=int(p1,x,0,u); c3=c2.*y1+c1.*(1-y1); y2=int(u.*p1,x,u,inf); y3=int(x.*p1,x,0,u); y=y2+y3; c4=c3./yc4=(10.*u/x+2200.+120.*u-120.*m)*(.50089433903710877308927024420163*erf(.35916569992135941929847473336317e-2*u-2.1549941995281565157908484001790)+.49973897687681912453571729041131)+(10.*u/x+1070.+4.*x)*(.50026102312318087546428270958869-.50089433903710877308927024420163*erf(.35916569992135941929847473336317e-2*u-2.1549941995281565157908484001790)/(.50089433903710877308927024420163*u-.50089433903710877308927024420163*erf(.35916569992135941929847473336317e-2*u-2.1549941995281565157908484001790)*u-78.682170542635657954656940854112*exp(-4.6440000000000000566656097045204-.12900000000000000157404471401446e-4*u2+.15480000000000000188885365681735e-1*u)+300.53660342226526385356214652098*erf(.35916569992135941929847473336317e-2*u-2.1549941995281565157908484001790)+300.60024125443500777024257422268)问题三的程序：syms x u m c;c1=0.98*0.98.*(2.*u./x.*10+1000)+(2*0.02*0.98+0.02*0.02).*(2.*u./x.*10+200.*x+3000+1500);c2=0.4*0.4.*(2.*u./x.*10+1000)+(2*0.4*0.6+0.6*0.6).*(2.*u./x.*10+3000+(u-m).*200);p1=0.00203*exp(-0.0000129.*(x-600).2);y1=int(p1,x,0,u);c3=c2.*y1+c1.*(1-y1);y2=int(u.*p1,x,u,inf);y3=int(x.*p1,x,0,u);y=y2+y3;c4=c3./yvpa(c4)c4=(20.*u./x+2680.+168.*u-168.*m).*(.50089433903710877308927024420163*erf(.35916569992135941929847473336317e-2.*u-2.1549941995281565157908484001790)+.49973897687681912453571729041131)+(20.*u./x+1138.6+7.92.*x).*(.50026102312318087546428270958869-.50089433903710877308927024420163.*erf(.35916569992135941929847473336317e-2.*u-2.1549941995281565157908484001790)./(.50089433903710877308927024420163.*u-.50089433903710877308927024420163.*erf(.35916569992135941929847473336317e-2.*u-2.1549941995281565157908484001790).*u-78.682170542635657954656940854112.*exp(-4.6440000000000000566656097045204-.12900000000000000157404471401446e-4.*u.2+.15480000000000000188885365681735e-1.*u)+300.53660342226526385356214652098.*erf(.35916569992135941929847473336317e-2.*u-2.1549941995281565157908484001790)+300.60024125443500777024257422268) x,u= meshgrid(20:40,340:360);m=0.95.*u;c4=(20.*u./x+2680.+168.*u-168.*m).*(.50089433903710877308927024420163*erf(.35916569992135941929847473336317e-2.*u-2.1549941995281565157908484001790)+.49973897687681912453571729041131)+(20.*u./x+1138.6+7.92.*x).*(.50026102312318087546428270958869-.50089433903710877308927024420163.*erf(.35916569992135941929847473336317e-2.*u-2.1549941995281565157908484001790)./(.50089433903710877308927024420163.*u-.50089433903710877308927024420163.*erf(.35916569992135941929847473336317e-2.*u-2.1549941995281565157908484001790).*u-78.682170542635657954656940854112.*exp(-4.6440000000000000566656097045204-.12900000000000000157404471401446e-4.*u.2+.15480000000000000188885365681735e-1.*u)+300.53660342226526385356214652098.*erf(.35916569992135941929847473336317e-2.*u-2.1549941995281565157908484001790)+300.60024125443500777024257422268)4.751750177348284.751947715119534.752772107047874.754168042472384.750250101955114.750407010970574.751190912601274.752546483969394.748855639394994.748972064645974.749715621912734.751030976017444.747566466005674.74764255086324.748345908083934.749621190106394.746382262309454.746418148543254.74708144843