cplex中文教程 第四章

合集下载

(完整word版)cplex翻译完全版

(完整word版)cplex翻译完全版

CPLEX 121. 简介................................................................. 3..2. 怎么用Cplex运行模型 (3)3. Cplex 概览........................................................... 3..3.1线性规划......................................................... 3.3.2二次约束规划....................................................4.3.3混合整数规划....................................................4.3.4可行松弛性...................................................... 5.3.5解池:产生和保持多解........................................... 5.4. GAMS 选项.......................................................... .9..5. Cplex选项总结 (10)5.1预处理和一般选项 (10)5.2单纯形法选项................................................... 1.25.3单纯形法的限制选项 (12)5.4单纯形法的容限选项 (13)5.5障碍特殊选项.................................................. 1.35.6筛选特殊选项................................................... 1.35.7混合整数规划选项 (13)5.8混合整数规划限制选项 (15)5.9混合整数规划解池选项 (16)5.10混合整数规划容许度选项....................................... 1.65.11输出选项...................................................... 1.75.12 GAMS/Cplex选项文件......................................... 1.76. 特殊备注 (18)6.1物理内存限制................................................... 1.86.2使用特殊有序集 (18)6.3使用半连续半整数变量 (19)6.4为求解MIP问题耗尽内存 (19)6.5不能证明整数最优 (20)6.6从混合整数规划的解开始 (20)6.7使用可行松弛性 (21)7. GAMS/ CPLE)日志文件 (22)8. CPLEX选项的详细说明 (25)1.简介GAMS/Cplex是一种用于GAMS (The General Algebraic Modeling System通用代数建模系统)的求解器,它使得用户可以把GAMS(通用代数建模系统的)的高级建模功能跟Cplex优化器的优势结合起来。

matlab-YAMLP-CPLEX调用与安装培训讲学

matlab-YAMLP-CPLEX调用与安装培训讲学
运行结果如图:
程序结果
算例2
目标函数: Max z=x1^2+x2^2+3*x3^2+4*x4^2+2*x5^2-8*x1-2*x2-3*x3-x4-2*x5 约束条件: 0<=xi<=99(i=1,2,...,5) x1+x2+x3+x4+x5<=400 x1+2*x2+2*x3+x4+6*x5<=800 2*x1+x2+6*x3<=800 x3+x4+5*x5<=200 在matlab中输入 x=intvar(1,5); f=[1 1 3 4 2]*(x'.^2)-[8 2 3 1 2]*x';F=set(0<=x<=99); F=F+set([1 1 1 1 1]*x'<=400)+set([1 2 2 1
1.定义变量:
sdpvar()表示实型; intvar()表示整型; binvar()表示0-1型;
2.确定目标函数f; 3.利用set设定约束条件(当有多个约束条件时,
可将各个约束条件用“+”连接) 4.求解目标函数solvesdp(F,-f)
最大值是-f,最小值是f。
5.利用double()查看相应结果。如double(f)等
6]*x'<=800)+set(2*x(1)+x(2)+6*x(3)<=800); F=F+set(x(3)+x(4)+5*x(5)<=200);solvesdp(F,-f) double(f) 80199 double(x) 53 99 99 99 0

CPLEX12.6在C++环境下的配置

CPLEX12.6在C++环境下的配置

CPLEX的安装
• 在安装CPLEX之前,应正确安装适用的编程平台 (例如Visual Studio、Java Eclipse)。 • 以CPLEX 12.61为例,点击安装程序,加载完成后 进行安装,基本步骤如下,一般选择默认路径(C 盘),中间步骤均选择默认即可。
CPLEX的安装
CPLEX的配置
求解简单线性规划
min x y z 2 x 3 y 10 4y+5z 15 s.t. 3x+z 11 x 0, 0 y 3, z 5
求解简单线性规划
求解人力资源调度问题
现有一个建造房子的项目,房子的建造一共包括5个主 要工作,分别为地基、墙壁、天花板、外墙装饰、内墙 装饰,且需要完成采购、监工、砖瓦工、油工、木工这 五种技能,每道工作各自需要其中的若干种技能,每种 所需技能都需要雇用多技能工人去完成。每道工作需要 占用一个固定的天数,为每道工作的执行时间。由于人 与人之间的差异性,并不是每个工人都有能力去执行每 项技能,他们都只掌握了执行其中部分技能的能力。并 且,由于实际情况的限制,各道工作之间存在着一定的 优先关系,例如,地基必须要先于其他各个工作之前完 成,墙壁必须要先于外墙装饰和内墙装饰之前完成等。 根据实际情况,本文给出如下限定: (1)房子的每个工作可能需要多个工人来进行,且这 几个工人必须同时开始,同时结束; (2)每个工人在每个单位时间至多只能执行一个技能 单位的工作,也就是说一个工人不能同时执行一项以上 的工作; (3)一项工作一旦开始,在完成之前不能有中断出现;
CPLEX的配置
CPLEX的配置
求解简单线性规划
min x y z 2 x 3 y 10 4y+5z 15 s.t. 3x+z 11 x 0, 0 y 3, z 5

Cplex学习笔记

Cplex学习笔记

Cplex学习笔记thedualsimplexmethodisthefirstchoiceforoptimizingalinearprogrammingproblem,esp eciallyforprimal-degenerateproblemswithlittlevariabilityintherighthandsidecoefficientsbutsignif icantvariabilityinthecostcoefficients.PrimarSimpleXOptimizer有时会更好地解决变量数量显著增加约束数量的问题,或者解决成本系数中几乎没有变化的问题。

六百五十四networkoptimizermayhaveapositiveimpactonperformanceofnetworkstructure 屏障优化因子和方法在大型稀疏问题上是有效的siftingwasdevelopedtoexploitthecharacteristicsofmodelswithlargeaspectratios(th atis,alargeratioofthenumberofcolumnstothenumberofrows).并行优化器多线程计算机平台简言之:原始变量多于约束对偶右停左变,边界大型稀疏阵筛选行列数差大。

要保存基本安全文件:libraryroutinecpxmbasewrite,afterthecalltotheoptimizer.然后从该文件中更新更新基础:cpxreadcopybase.如果解决问题所需的迭代次数与不足的次数大致相同,那么你做得很好。

如果迭代次数比行数(或更多)大三倍,则完全有可能通过改变为UALSimplexOptimizer或Primind implexOptimizer设置定价算法的参数DPRI来提高性能ifyouobservethatyourproblemhasdifficultystayingfeasibleduringitssolution,theny oushouldconsideranalternativescalingmethod.scaindparametersettingsforscalingme thods斯坎德价值意义-1noscaling0平衡缩放(默认)1缩放缩放671cplex学习手册STL标准模板库简称STL,即标准模板库。

CPLEX中文教程(第七章)

CPLEX中文教程(第七章)

CPLEX中文教程(第七章)第七章:IBMILOGCPLE某在公交乘务排班优化问题中的应用(一)公交乘务排班优化问题概述公交乘务排班问题属活动资源的优化利用问题。

一般是根据给定的乘务任务、乘务规则等条件,考虑一定的优化目标,对乘务员(组)的出乘时间、地点,担当的乘务任务、时刻,退乘时间、地点等做出具体安排,以确保一定周期内的所有乘务任务被执行。

第七章:IBMILOGCPLE某在公交乘务排班优化问题中的应用公交乘务排班问题的基本元素和各元素的基本属性如下:1)公交线路:具有出发站、出发时刻、到达站、到达时刻、中途停站等基本属性;2)乘务员类型(组):包括司机、售票员等属性;3)乘务规则:包括间休时间、工作时间、休息时间、乘务周期、月工时等乘务值乘规则。

4)目标函数:乘务成本最小、需要的乘务员数量最少等。

5)约束条件:乘务员的工作时间必须满足乘务规则;每个线路、每个班次都必须有乘务员值乘;乘务员劳动负衡均衡等约束条件。

第七章:IBMILOGCPLE某在公交乘务排班优化问题中的应用例1.下面是某条线路的基本情况:1、该线路的开收班时间:冬令(12月~3月):6:20~18:10,夏令(4月~11月):6:15~18:202、该线路的司机人数:15人3、该线路排班间隔:平时:8~10分钟/班;上下班高峰(6:00~8:30,11:30~13:30,16:30~18:00):4~8分钟/班节假日:5~10分钟/班4、该线路的运行时间:正常:80~85分钟/班高峰:100~120分钟/班规定:(1)司机每天上班时间不超过8小时;(2)司机连续开车不得超过4小时;(3)每名司机至少每月完成120班次。

问题一:针对五月份的节假日和非节假日,分别求出每日最少班次总数;问题二:阐述你对上述规定的理解,并根据你的理解建立适当的数学模型,合理地设计五月份该线路的司机排班方案。

第七章:IBMILOGCPLE某在公交乘务排班优化问题中的应用(二)模型建立与求解1、问题一:数学模型的建立(1)假设将非节假日一天的工作时间分为n个时段6minZ=i=1某ikk某imiQii1i1kkS.T.某imimkQii1i18mi104mi8k1,2,3,4,5,6.k1,2,3,4,5,6.i2,4,6.i1,3,5.Qi——第i个时段的时长,i=1,2…,nmi——表示发车间隔;某i——表示发车班次用第七章:IBMILOGCPLE某在公交乘务排班优化问题中的应用问题一:1、数学模型的建立(2)若节假日,用Q表示节假日一天的工作时长,发车间隔用m表示,则例1中节假日最少班次的数学公式为:QminZ=,m5≤m≤10.第七章:IBMILOGCPLE某在公交乘务排班优第七章:IBMILOGCPLE某在公交乘务排班优化问题中的应用问题一的模型文件(续):某["T1"]某m["T1"]+某["T2"]某m["T2"]+(某["T3"]-1)某m["T3"]<=Q["T1"]+Q["T2"]+Q["T3"];Q["T1"]+Q["T2"]+Q["T3"]<=某["T1"]某m["T1"]+某["T2"]某m["T2"]+某["T3"]某m["T3"];某["T1"]某m["T1"]+某["T2"]某m["T2"]+某["T3"]某m["T3"]+(某["T4"]1)某m["T4"]<=Q["T1"]+Q["T2"]+Q["T3"]+Q["T4"];Q["T1"]+Q["T2"]+Q["T3"] +Q["T4"]<=某["T1"]某m["T1"]+某["T2"]某m["T2"]+某["T3"]某m["T3"]+某["T4"]某m["T4"];某["T1"]某m["T1"]+某["T2"]某m["T2"]+某["T3"]某m["T3"]+某["T4"]某m["T4"]+(某["T5"]1)某m["T5"]<=Q["T1"]+Q["T2"]+Q["T3"]+Q["T4"]+Q["T5"];Q["T1"]+Q["T2"] +Q["T3"]+Q["T4"]+Q["T5"]<=某["T1"]某m["T1"]+某["T2"]某m["T2"]+某["T3"]某m["T3"]+某["T4"]某m["T4"]+某["T5"]某m["T5"];某["T1"]某m["T1"]+某["T2"]某m["T2"]+某["T3"]某m["T3"]+某["T4"]某m["T4"]+某["T5"]某m["T5"]+(某["T6"]1)某m["T6"]<=Q["T1"]+Q["T2"]+Q["T3"]+Q["T4"]+Q["T5"]+Q["T6"];Q["T1"] +Q["T2"]+Q["T3"]+Q["T4"]+Q["T5"]+Q["T6"]<=某["T1"]某m["T1"]+某["T2"]某m["T2"]+某["T3"]某m["T3"]+某["T4"]某m["T4"]+某["T5"]某m["T5"]+某["T6"]某m["T6"];}第七章:IBMILOGCPLE某在公交乘务排班优化问题中的应用求解结果:高峰期(6:158:30)发车时间班次1234567891011121314151617平时(8:3011:30)发车时间08:3108:4108:5109:0109:1109:2109:3109:4109:5110:0110:1110:2110: 3110:4110:5111:0111:1111:21班次181920222223242526272829303132333435高峰期(11:3013:30)发车时间11:3111:3911:4711:5512:0312:1112:1912:2712:3512:4312:5112:5913: 0713:1513:23班次363738394041424344454647484950平时(13:3016:30)发车时间13:3113:4113:5114:0114:1114:2114:3114:4114:5115:0115:1115:2115: 3115:4115:5116:0116:1116:21班次515253545556575859606162636465666768高峰期(16:3018:00)发车时间16:3116:3916:4716:5517:0317:1117:1917:2717:3517:4317:5117:59班次697071727374757677787980平时(18:0018:20)发车时间18:0718:17班次8182//olution(optimal)withobjective82某=[17181518122]06:1506:2306:3106:3906:4706:5507:0307:1107:19据此,得到该线路非节假日的最少班次排班如右:07:2707:3507:4307:5107:5908:0708:1508:23第七章:IBMILOGCPLE某在公交乘务排班优化问题中的应用求解结果:发车时间班次123456789101112发车时间08:1508:2508:3508:4508:5509:0509:1509:2509:3509:4509:5510:05班次131415161718192022222324时间(6:20-18:10)发车时间10:1510:2510:3510:4510:5511:0511:1511:2511:3511:4511:5512:05班次252627282930313233343536发车时间12:1512:2512:3512:4512:5513:0513:1513:2513:3513:4513:5514:05班次373839404142434445464748发车时间14:1514:2514:3514:4514:5515:0515:1515:2515:3515:4515:5516:05班次495051525354555657585960发车时间16:1516:2516:3516:4516:5517:0517:1517:2517:3517:4517:5518:0518:15班次61626364656667686970717273节假日的最少班次排班如右:06:1506:2506:3506:4506:5507:0507:1507:2507:3507:4507:5508:05第七章:IBMILOGCPLE某在公交乘务排班优化问题中的应用问题二:1、数学模型的建立(1)非节假日目标函数:min约束条件:某jj119式中:ai——表示非节假日一天内第i个司机工作的高峰班次数bi——表示非节假日一天内第i个司机工作的非高峰班次数某j——表示司机排班情况(aij119j119某j)44(bi某j)38第七章:IBMILOGCPLE某在公交乘务排班优化问题中的应用司机排班情况第七章:IBMILOGCPLE某在公交乘务排班优化问题中的应用模型编码:{tring}Categorie={"C1","C2","C3","C4","C5","C6","C7","C8","C9","C10","C11","C12","C13","C14","C15","C16","C17","C18","C19"};inta [Categorie]=[0,0,0,0,0,1,1,1,1,1,2,2,2,2,3,3,3,4,4];intb[Categor ie]=[1,2,3,4,5,0,1,2,3,4,0,1,2,3,0,1,2,0,1];dvarint+某[Categorie]in0..50;minimizeum(oinCategorie)某[o];ubjectto{um(oinCategorie)a[o]某某[o]==44;um(oinCategorie)b[o]某某[o]==38;}第七章:IBMILOGCPLE某在公交乘务排班优化问题中的应用求解结果://olution(optimal)withobjective17某=[01005000000000000011]该结果给出了一天需要的最少司机人数为17人,同时给出了一个可行的排班方案,即17位司机中,有1位采用第2种排班情况,有5位采用第5种排班情况,有11位采用第19种排班情况。

cplex中文教程 第五章

cplex中文教程 第五章
A4 A3 A2 A1 A0
T集 t节
600 N40=60 N41=250 N42=180
550 2.5
500 3.0 2.0
N30=130 N31=100 N20=300
5个支点方向的直达车流图
第五章 IBM ILOG CPLEX在车流组织优化中的应用 (二) 货物列车编组计划的整数规划模型 变量定义: ������ (1)������������������ ——表示车流(i, j)(含������������ 后方站发往������������ 站改 ������ 编中转的车流)在������������ 站改编的车流量,这里i>j>k,������������������ 为非负整数; (2)������������������ ——0-1变量,表示是否开行列流,若开行,则 ������������������ = 1,反之������������������ = 0,这里i>j。 (3)������ ——所有直达列流在始发站产生的集结车小时 集 总消耗; (4)������ ——所有直达车流在图中支点站的改编车小时 改 总消耗;
第五章 IBM ILOG CPLEX在车流组织优化中的应用 求解结果:
即:min������ = 2360车 · h,根据运行结果可得出最优货物列车 耗 编组计划图:
A4
A3
A2
A1
A0
TБайду номын сангаас t节
600
550 2.5 60+250
500 3.0 2.0
130+300 N43+180 N32+180+130+130 N21+100 N10+60

【CPLEX教程02】配置Cplex的Java环境以及API说明

【CPLEX教程02】配置Cplex的Java环境以及API说明

【CPLEX教程02】配置Cplex的Java环境以及API说明因为⼩编⼀般⽤的C++和Java⽐较多,⽽且现在开发⼤型算法⽤这类⾯向对象的编程语⾔也⽅便得多。

基于上⾯的种种考虑,加上时间和精⼒有限,所以就暂时只做C++和Java的详细教程辣。

关于matlab和python的也许后续会补上的吧。

然后在开始之前,照例先把环境给配置好。

那么就先配置java的环境吧。

CPLEX系列教程可以关注我们的公众号哦!获取更多精彩消息!前⾯已经说了怎么下载和安装cplex了,如图:确保已经安装上这个版本,我们才能开始下⼀步的⼯作。

java⼩编⼀般⽤的ide是eclipse,就配置⼀下关于eclipse的。

其他的开发环境请⼤家⾃⾏设置哈。

新建⼀个⼯程,添加⼀个package,添加⼀个带main函数的类。

代码先别写。

在项⽬右键,选择build path -> Configure Build Path……找到Libraries->Add External JARs……,然后定位到\lib这个⽂件夹,把cplex.jar给添加进去。

到这⼀步还不⾏,还需要把CPLEX的动态运⾏库给添加进去,好让java程序运⾏的时候能够找到,具体做法是:在项⽬右键,选择build path -> Configure Build Path……,找到Libraries,点开JRE System Library,在Native library location那⾥点edit,把cplex下的\bin\x64_win64⽂件夹给添加进去,这⾥⾯有程序运⾏所需要动态库。

⼀个简单的线性规划问题:把下⾯代码复制进main函数⾥⾯:try {IloCplex cplex = new IloCplex(); // creat a modeldouble[] lb = {0.0, 0.0, 0.0};double[] ub = {40.0, Double.MAX_VALUE, Double.MAX_VALUE};IloNumVar[] x = cplex.numVarArray(3, lb, ub);double[] objvals = {1.0, 2.0, 3.0};cplex.addMaximize(cplex.scalProd(x, objvals));double[] coeff1 = {-1.0, 1.0, 1.0};double[] coeff2 = {1.0, -3.0, 1.0};cplex.addLe(cplex.scalProd(x, coeff1), 20.0);cplex.addLe(cplex.scalProd(x, coeff2), 30.0);if (cplex.solve()) {cplex.output().println("Solution status = " + cplex.getStatus());cplex.output().println("Solution value = " + cplex.getObjValue());double[] val = cplex.getValues(x);for (int j = 0; j < val.length; j++)cplex.output().println("x" + (j+1) + " = " + val[j]);}cplex.end();} catch (IloException e) {System.err.println("Concert exception caught: " + e);}点击运⾏,出现下⾯的结果就成功啦:最后,如果提⽰找不到build path ,share libraries什么的,请确保第⼆步配置正确!⾄此,我们已经能愉快使⽤cplex啦。

第4章-方程求解(Maple中文教程)

第4章-方程求解(Maple中文教程)

第4章-⽅程求解(Maple中⽂教程)第四章⽅程求解1 代数⽅程(组)求解1.1 常⽤求解⼯具—solve求解代数⽅程或代数⽅程组, 使⽤Maple 中的solve 函数. 求解关于x 的⽅程eqn=0的命令格式为:solve(eqn, x);求解关于变量组vars 的⽅程组eqns 的命令为:solve(eqns, vars);> eqn:=(x^2+x+2)*(x-1);:= eqn () + + x 2x 2() ? x 1> solve(eqn,x);,,1? + 1212I 7? ? 1212I 7 当然, solve 也可以求解含有未知参数的⽅程:> eqn:=2*x^2-5*a*x=1;:= eqn = ? 2x 25a x 1 > solve(eqn,x);, + 54a 14 + 25a 28 ? 54a 14 + 25a 28 solve 函数的第⼀个参数是有待求解的⽅程或⽅程的集合, 当然也可以是单个表达式或者表达式的集合, 如下例:> solve(a+ln(x-3)-ln(x),x);3e a+ 1ea 对于第⼆个参数, Maple 的标准形式是未知变量或者变量集合, 当其被省略时, 函数indets ⾃动获取未知变量. 但当⽅程中含有参数时, 则会出现⼀些意想不到的情况: > solve(a+ln(x-3)-ln(x));{}, = x x = a ? + ()ln ? x 3()ln x很多情况下, 我们知道⼀类⽅程或⽅程组有解, 但却没有解决这类⽅程的⼀般解法, 或者说没有解析解. ⽐如, ⼀般的五次或五次以上的多项式, 其解不能写成解析表达式. Maple 具备⽤所有⼀般算法尝试所遇到的问题, 在找不到解的时候, Maple 会⽤RootOf 给出形式解.> x^7-2*x^6-4*x^5-x^3+x^2+6*x+4;+ + + x 72x 64x 5x 3x 26x 4> solve(%); + 15 ? 15()RootOf , ? ? _Z 5_Z 1 = index 1()RootOf , ? ? _Z 5_Z 1 = index 2(RootOf ,) ? ? _Z 5_Z 1 = index 3,,,, ()RootOf , ? ? _Z 5_Z 1 = index 4()RootOf , ? ? _Z 5_Z 1 = index 5,,> solve(cos(x)=x,x);()RootOf ? _Z ()cos _Z对于⽅程组解的个数可⽤nops 命令获得, 如:> eqns:={seq(x[i]^2=x[i],i=1..7)};:= eqns {,,,,,, = x 12x 1 = x 22x 2 = x 32x 3 = x 42x 4 = x 52x 5 = x 62x 6 = x 72x 7} > nops({solve(eqns)});128但是, 有时候, Maple 甚⾄对⼀些“显⽽易见”的结果置之不理, 如:> solve(sin(x)=3*x/Pi,x);()RootOf ? 3_Z ()sin _Z π此⽅程的解为0 ,6π±, 但Maple 却对这个超越⽅程⽆能为⼒, 即便使⽤allvalues求解也只有下述结果:> allvalues(%);()RootOf , ? 3_Z ()sin _Z π0.另外⼀个问题是, Maple 在求解⽅程之前,会对所有的⽅程或表达式进⾏化简, ⽽不管表达式的类型, 由此⽽产⽣⼀些低级的错误: > (x-1)^2/(x^2-1);() ? x 12x 21> solve(%);1但是, ⼤量实验表明, solve 的确是⼀个实⽤的⽅程求解⼯具, 但是也不可盲⽬相信它给出的⼀切结果, 特别是对于⾮线性⽅程⽽⾔, 对于给出的结果需要加以验证.下⾯通过⼏个例⼦说明在Maple 中⾮线性⽅程组的求解问题.例:求解⽅程组: ?=?=+y x y x 925222> eqns:={x^2+y^2=25,y=x^2-5};:= eqns {}, = y ? x 25 = + x 2y 225> vars:={x,y};:= vars {},x y> solve(eqns,vars);,,,{}, = x 0 = y -5{}, = x 0 = y -5{}, = y 4 = x 3{}, = y 4 = x -3也可⽤下⾯的语句⼀步求出:> solve({x^2+y^2=25,y=x^2-5},{x,y});,,,{}, = x 0 = y -5{}, = x 0 = y -5{}, = y 4 = x 3{}, = y 4 = x -3这个问题⾮常简单, 但通常遇到的⾮线性问题却不是这么简单, 例如要求解⽅程组:y x y x y x ?=+=+,122> eqns:={x^2+y^2=1,sqrt(x+y)=x-y}; vars:={x,y};:= eqns {}, = + x 2y 21 = + x y ? x y:= vars {},x y> sols:=solve(eqns,vars);sols = y ()RootOf , + + 2_Z 24_Z 3 ? -1.000000000.7071067812I ,{ := = x ? ? ()RootOf , + + 2_Z 24_Z 3 ? -1.000000000.7071067812I 2}{}, = x 1 = y 0,可以看出, ⽅程解的形式是以集合的序列给出的, 序列中的每⼀个集合是⽅程的⼀组解, 这样就很利于我们⽤subs 把解代⼊原⽅程组进⾏检验:> subs(sols[2],eqns);{} = 11> sols2:=allvalues(sols[1]);:= sols2{}, = x ? + 11I 2 = y ? ? 11I 2 > simplify(subs(sols2,eqns));{}, = I 2I 2 = 111.2 其他求解⼯具1.2.1 数值求解对于求代数⽅程的数值解问题, Maple 提供了函数fsolve , fsolve 的使⽤⽅法和solve 很相似:fsolve(eqns, vars, options);其中, eqns 表⽰⼀个⽅程、⽅程组或者⼀个程序, vars 表⽰⼀个未知量或者未知量集合, options 控制解的参数(诸如:complex:复根; maxsols=n :只找到n 阶最⼩根; intervals :在给定闭区间内求根, 等).> fsolve(x^5-x+1,x);-1.167303978> fsolve(x^5-x+1,x,complex);-1.167303978 ? -.1812324445 1.083954101I + -.1812324445 1.083954101I ? .7648844336.3524715460I ,,, +.7648844336.3524715460I ,> fsolve(x^3-3*x+1,x,0..1);.3472963553对于多项式⽅程, fsolve 在默认情况下以给出所有的实数解, 如果附加参数complex , 就可以给出所有的解. 但对于更⼀般的其他形式的⽅程, fsolve 却往往只满⾜于得到⼀个解:> eqn:=sin(x)=x/2;:= eqn = ()sin x 12x > fsolve(eqn);0.> fsolve(eqn,x,0.1..infinity);1.895494267> fsolve(eqn,x,-infinity..-0.1);-1.895494267函数fsolve 主要基于两个算法, 通常使⽤⽜顿法, 如果⽜顿法⽆效, 它就改⽽使⽤切线法. 为了使fsolve 可以求得所有的实根, 我们通常需要确定这些根所在的区间. 对于单变量多项式, 函数realroot 可以获得多项式的所有实根所在的区间.> 4+6*x+x^2-x^3-4*x^5-2*x^6+x^7;+ + ? ? ? + 46x x 2x 34x 52x 6x 7> realroot(%);[],,[],02[],24[],-2-1函数realroot 还有⼀个可选参数, 它是⽤来限制区间的最⼤长度的, 为了保证使⽤数值求解⽅法时收敛, 我们可以⽤它限制区间的最⼤长度:> realroot(%%,1/1000);,,,1195299,33131657,-633-1265 求解⽅程或⽅程组的整数解时使⽤函数isolve , 它常常被⽤来求解不定⽅程. 例如著名的“百钱买百鸡”问题?的求解过程为:> isolve({x+y+z=100,5*x+3*y+z/3=100});{},, = z + 753_Z1 = x 4_Z1 = y ? 257_Z1据此可得满⾜该问题的三组解为:{x, y, z}={4, 18, 78}, {x, y, z}={8, 11, 81}, {x, y, z}={12, 4, 84}1.2.2 整数环中的⽅程(组)求解利⽤Maple 中的函数msolve(eqns, vars, n), 可以在模n 的整数环中求解⽅程(组)eqns.例:在Z 7中求解Pell ⽅程2837?=x y > msolve(y^7=x^3-28,7);{}, = x 3 = y 6{}, = x 4 = y 1{}, = y 0 = x 0{}, = x 1 = y 1{}, = y 6 = x 6,,,,,{}, = x 2 = y 1{}, = y 6 = x 5, 再如下例:> msolve(y^4=x^3+32,5);,,,,{}, = x 2 = y 0{}, = x 4 = y 1{}, = x 4 = y 2{}, = x 4 = y 3{}, = x 4 = y 41.2.3 递归⽅程的求解在Maple 中, 可以求解有限差分⽅程(也称递归⽅程), 所需调⽤的函数是rsolve , 该函数使⽤的是⼀些⽐较通⽤的⽅法, 例如产⽣函数法、z 变换法以及⼀些基于变量替换和特征⽅程的⽅法. 作为例⼦, 求解Fibonacci 多项式:> eq:=f(n)=f(n-1)+2*f(n-2);:= eq = ()f n + ()f ? n 12()f ? n 2> rsolve({eq,f(0)=1,f(1)=1},f(n)); + 13()-1n 232n 当然, 并不是所有的递归形式的函数⽅程的解可以写成解析形式, 如果不能, Maple 将保留原来的调⽤形式. 此时, 可⽤asympt 函数获得它的渐进表达式, 也就是1/n 的级数解. 例如, 对于⼀个具有超越形式的递归函数⽅程, 仍然可以得到解的渐进形式: ? 百钱买百鸡问题:⽤100元钱买100只鸡, ⼤公鸡5元钱1只, ⼤母鸡3元钱1只, ⼩鸡1元钱3只, 问如何买法?> rsolve(u(n+1)=ln(u(n)+1),u(n));()rsolve , = ()u + n 1()ln + ()u n 1()u n> asympt(%,n,5);+ + + 21n + _C 23()ln n n 2 ? + ? + 1913_C 12_C 2 + 23_C 29()ln n 29()ln n 2n 3????????O 1n 41.2.4 不等式(组)求解求解⼀元不等式⽅程(组)使⽤命令solve :> solve((x-1)*(x-2)*(x-3)<0,x);,()RealRange ,?∞()Open 1()RealRange ,()Open 2()Open 3> solve((x-1+a)*(x-2+a)*(x-3+a) < 0, {x});,{} < x ? 1a {}, < ? 2a x < x ? 3a> solve(exp(x)>x+1);,()RealRange ,?∞()Open 0()RealRange ,()Open 0∞> solve({x^2*y^2=0,x-y=1,x<>0});,{}, = y 0 = x 1{}, = y 0 = x 1对于由不等式⽅程组约束的最优问题的求解使⽤“线性规则”⼯具包simplex : > with(simplex):> cnsts:={3*x+4*y-3*z<=23, 5*x-4*y-3*z<=10,7*x+4*y+11*z<=30};:= cnsts {,, ≤ } + ? 3x 4y 3z 23 ≤ ? ? 5x 4y 3z 10 ≤ + + 7x 4y 11z 30 > obj:=-x+y+2*z;:= obj ? + + x y 2z> maximize(obj,cnsts union {x>=0,y>=0,z>=0});{},, = z 12 = y 498= x 0 2 常微分⽅程求解微分⽅程求解是数学研究与应⽤的⼀个重点和难点. Maple 能够显式或隐式地解析地求解许多微分⽅程求解. 在常微分⽅程求解器dsolve 中使⽤了⼀些传统的技术例如laplace 变换和积分因⼦法等, 函数pdesolve 则使⽤诸如特征根法等经典⽅法求解偏微分⽅程. 此外, Maple 还提供了可作摄动解的所有⼯具, 例如Poincare-Lindstedt 法和⾼阶多重尺度法.帮助处理常微分⽅程(组)的各类函数存于Detools 软件包中, 函数种类主要有:可视化类的函数, 处理宠加莱动态系统的函数, 调整微分⽅程的函数, 处理积分因⼦、李对称法和常微分⽅程分类的函数, 微分算⼦的函数, 利⽤可积性与微分消去的⽅法简化微分⽅程的函数, 以及构造封闭解的函数等. 更重要的是其提供的强⼤的图形绘制命令Deplot 能够帮助我们解决⼀些较为复杂的问题.2.1 常微分⽅程的解析解求解常微分⽅程最简单的⽅法是利⽤求解函数dsolve . 命令格式为:dsolve(ODE);dsolve(ODE, y(x), extra_args);dsolve({ODE, ICs}, y(x), extra_args);dsolve({sysODE, ICs}, {funcs}, extra_args);其中, ODE—常微分⽅程, y(x)—单变量的任意变量函数, Ics—初始条件, {sysODE}—ODE ⽅程组的集合, {funcs}—变量函数的集合, extra_args—依赖于要求解的问题类型.例如, 对于⼀阶常微分⽅程y xy y y x ?=′)ln(可⽤dsolve 直接求得解析解: > ODE:=x*diff(y(x),x)=y(x)*ln(x*y(x))-y(x);:= ODE = xx ()y x ? ()y x ()ln x ()y x ()y x > dsolve(ODE,y(x));= ()y x ex可以看出, dsolve 的第⼀个参数是待求的微分⽅程, 第⼆个参数是未知函数. 需要注意的是, ⽆论在⽅程中还是作为第⼆个参数,未知函数必须⽤函数的形式给出(即:必须加括号, 并在其中明确⾃变量), 这⼀规定是必须的, 否则Maple 将⽆法区分⽅程中的函数、⾃变量和参变量, 这⼀点和我们平时的书写习惯不⼀致. 为了使其与我们的习惯⼀致, 可⽤alias 将函数⽤别称表⽰:> alias(y=y(x));> ODE:=x*diff(y,x)=y*ln(x*y)-y;:= ODE = xx y ? y ()ln x y y > dsolve(ODE,y);= y ex _C1x函数dsolve 给出的是微分⽅程的通解, 其中的任意常数是⽤下划线起始的内部变量表⽰的.在Maple 中, 微分⽅程的解是很容易验证的, 只需要将解代⼊到原⽅程并化简就可以了.> subs(%,ODE);= xx e ??x _C1x ? e ??x _C1??ln e ??x _C1x e ??x_C1x > assume(x,real): assume(_C1,real):> simplify(%); = ?e x~_C1~()? + x~_C1~x~_C1~?e x~_C1~()? + x~_C1~x~_C1~> evalb(%);trueevalb 函数的⽬的是对⼀个包含关系型运算符的表达式使⽤三值逻辑系统求值, 返回的值是true, false 和FAIL. 如果⽆法求值, 则返回⼀个未求值的表达式. 通常包含关系型运算符“=, <>, <, <=, >, >=”的表达式在Maple 中看作是代数⽅程或者不等式. 然⽽, 作为参数传递给evalb 或者出现在if 或while 语句的逻辑表达式中时, 它们会被求值为true 或false. 值得注意的是, evalb 不化简表达式, 因此在使⽤evalb 之前应将表达式化简, 否则可能会出错. 再看下⾯常微分⽅程的求解:12+=′y y> alias(y=y(x)):> ODE:=diff(y,x)=sqrt(y^2+1);:=ODE = ??xy + y 21 > dsolve(ODE,y); = y ()sinh + x _C1函数dsolve 对于求解含有未知参变量的常微分⽅程也完全可以胜任:> alias(y=y(x)):> ODE:=diff(y,x)=-y/sqrt(a^2-y^2);:=ODE = ?y ?y ? a 2y 2> sol:=dsolve(ODE,y); := sol = + ? + x ? a 2y 2a 2ln + 2a 22a 2a 2y 2y a 2_C10由此可见, 对于不能表⽰成显式结果的微分⽅程解, Maple 尽可能将结果表⽰成隐式解. 另外, 对于平凡解y=0常常忽略, 这⼀点应该引起注意.dsolve 对于求解微分⽅程初值问题也⼗分⽅便的:> ODE:=diff(u(t),t$2)+omega^2*u(t)=0;:= ODE = +2t 2()u t ω2()u t 0 > dsolve({ODE,u(0)=u0,D(u)(0)=v0},u(t)); = ()u t + v0()sin ωt ωu0()cos ωt 2.2 利⽤积分变换求解微分⽅程对于特殊的微分⽅程, 我们还可以指定dsolve 利⽤积分变换⽅法求解, 只需要在dsolve 中加⼊可选参数method=transform 即可.其中transform 是积分变换, 可以是laplace 、fourier 、fouriercos 或者fouriersin 变换.作为例⼦, 我们来看⼀个具有阻尼的振⼦在阶跃冲击(Heaviside 函数)下的响应: >ODE:=diff(u(t),t$2)+2*d*omega*diff(u(t),t)+omega^2*u(t)=Heaviside(t);:= ODE = + + 2t 2()u t 2d ωt ()u t ω2()u t ()Heaviside t > initvals:=(u(0)=u[0],D(u)(0)=v[0]);:= initvals , = ()u 0u 0 = ()()D u 0v 0> solution:=dsolve({ODE,initvals},u(t),method=laplace);:= solution = ()u t + 1ωe ()?t d ω + () ? ω2u 01()cosh t ? d 2ω2ω2ω() + ? ωv 0d ω2u 0d ()sinh t ? d 2ω2ω2 ? d 2ω2ω2ωMaple 给出了问题的通解, 但没有区分⾃由振动(d=0)、⽋阻尼(01)的情况. 下⾯加以区分求解:> assume(omega>0):> simplify(subs(d=0,solution));= ()u t + ? + 1()cos t ωω2u 0()cos t ωv 0()sin t ωωω2> K:=subs(d=1/5,u[0]=1,v[0]=1,solution);:= K = ()u t + 1ωe ()?/15t ω + () ? ω21cosh t ?2425ω2ω + ? ω15ω215 sinh t ?2425ω2?2425ω2> with(plots):> plot3d(rhs(%%),omega=2/3..4/3,t=0..20,style=hidden,orientation=[-30,45],axes=framed);对于d=1的情况, 可可⽤下式获得结果:> limit(rhs(solution),d=1);() + ? + ? + ω2u 0ω2v 0t 1ω3u 0t t ωe()t ωe ()?t ωω2再如下例:> diff(u(t),t$2)+3*diff(u(t),t)+2*u(t)=exp(-abs(t));= + +2t 2()u t 3?()u t 2()u t e ()?t > dsolve(%,u(t),method=fourier);= ()u t + + ? 23e ()?2t ()Heaviside t 16e t ()Heaviside ?t e ()?t t ()Heaviside t 12e ()?t ()Heaviside t2.3 常微分⽅程组的求解函数dsolve 不仅可以⽤来求解单个常微分⽅程, 也可以求解联⽴的常微分⽅程组. 特别是对于线性微分⽅程组, 由于数学上具有成熟的理论, Maple 的求解也是得⼼应⼿. 其命令格式为:dsolve( {eqn1, eqn2, …, ini_conds}, {vars});其中, ini_conds 是初始条件.> eqn1:={diff(x(t),t)=x(t)+y(t),diff(y(t),t)=y(t)-x(t)};:= eqn1{, = }??t ()x t + ()x t ()y t = ??t()y t ? ()y t ()x t > dsolve(eqn1,{x(t),y(t)});{, = ()x t e }t () + _C1()sin t _C2()cos t = ()y t e t () ? _C1()cos t _C2()sin t >eqn2:=2*diff(x(t),t$2)+2*x(t)+y(t)=2*t;:= eqn2 = + + 22t 2()x t 2()x t ()y t 2t > eqn3:=diff(y(t),t$2)+2*x(t)+y(t)=t^2+1;:= eqn3 = + +2t 2()y t 2()x t ()y t + t 21 > dsolve({eqn2, qn3, x(0)=0, D(x)(0)=1, y(0)=0, D(y)(0)=0}, {x(t),y(t)} ); = ()x t + ? + 18()sin 2t 2112t 3148t 434t ,{ = ()y t ? + ? + 14()sin 2t 212t 12t 216t 3124t 4} 2.4 常微分⽅程的级数解法1) 泰勒级数解法当⼀个常微分⽅程的解析解难以求得时, 可以⽤Maple 求得⽅程解的级数近似, 这在⼤多数情况下是⼀种⾮常好的⽅法. 级数解法是⼀种半解析半数值的⽅法. 泰勒级数法的使⽤命令为:dsolve({ODE,Ics}, y(x), 'series'); 或dsolve({ODE,Ics}, y(x), 'type=series');下⾯求解物理摆的⼤幅振动⽅程:θθsin g l ?=, 其中l 是摆长,θ是摆⾓, g 是重⼒加速度.> ODE:=l*diff(theta(t),t$2)=-g*sin(theta(t)); := ODE = l2t 2()θt ?g ()sin ()θt > initvals:=theta(0)=0,D(theta)(0)=v[0]/l; := initvals , = ()θ00 = ()()D θ0v 0l> sol:=dsolve({ODE,initvals},theta(t),type=series);:= sol = ()θt ? + + v 0l t 16g v 0l 2t 31120g v 0() + v 02g l l4t 5()O t 6 > Order:=11:> sol:=dsolve({ODE,initvals},theta(t),type=series); := sol = ()θt ? + ? + + v 0t 1g v 0l 2t 31g v 0() + v 02g l l 4t 51g v 0() + + 11g l v 02g 2l 2v 04l 6t 71g v 0() + + + 57g v 04l 102g 2v 02l 2g 3l 3v 06l 8t9()O t 112) 幂级数解法对于⼀个符号代数系统来说, 幂级数是必不可少的微分⽅程求解⼯具. 幂级数求解函数powsolve 存于⼯具包powseries 中. 但是, 这⼀求解函数的使⽤范围很有限, 它只可以⽤来求解多项式系数的线性常微分⽅程或⽅程组,其求解命令为:powseries[function] (prep)或直接载⼊软件包后⽤function(prep), prep 为求解的线性微分⽅程及其初值.例:求解:042=+′′+′y x y y x > ODE:=x*diff(y(x),x$2)+diff(y(x),x)+4*x^2*y(x)=0;:= ODE = + + x2x 2()y x ??x ()y x 4x 2()y x 0 > dsolve(ODE,y(x));= ()y x + _C1BesselJ ,043x ()/32_C2??????BesselY ,043x ()/32 > initvals:=y(0)=y0,D(y)(0)=0;:= initvals , = ()y 0y0 = ()()D y 00> with(powseries):> sol:=powsolve({ODE,initvals});:= sol proc ()... end proc powparm> tpsform(sol,x,16);+ + + y049y0x 3481y0x 6166561y0x 9459049y0x 121613286025y0x 15()O x 16 也可以⽤powsolve 给出的函数直接获得⽤递归形式定义的幂级数系数, 不过参数必须⽤_k , 这是powsolve 使⽤的临时变量.> sol(_k);4()a _k 3_k2 例:求解⼀维谐振⼦的解:0)(2=?+′′y x y ε> alias(y=y(x)):> ODE:=diff(y,x$2)+(epsilon-x^2)*y=0;:= ODE = +2x 2y () ? εx 2y 0 > H:=powsolve(ODE);:= H proc ()... end proc powparm> tpsform(H,x,8);C0C1x 12εC0x 216εC1x 3 + 124ε2C0112C0x 4+ 1120ε2C1120C1x 5????????? ? 130ε + 124ε2C0112C0160εC0 + ? ? + + + x 6????????? ? 1ε + 1ε2C11C11εC1x 7()O x 8 + +> H(_k);ε()a ? _k 2(a ) ? _k 4_k ()_k 1 2.5 常微分⽅程的数值解法在对微分⽅程的解析解失效后, 可以求助于数值⽅法求解微分⽅程. 数值求解的好处是只要微分⽅程的条件⾜够多时⼀般都可求得结果, 然⽽所得结果是否正确则必须依赖相关数学基础加以判断. 调⽤函数dsolve 求常微分⽅程初值问题的数值解时需加⼊参数type=numeric .另⼀⽅⾯, 常微分⽅程初值问题数值求解还可以选择算法, 加⼊参数“method=⽅法参数”即可, ⽅法参数主要有:rkf45:4~5阶变步长Runge-Kutta-Fehlberg 法dverk78:7~8阶变步长Runge-Kutta-Fehlberg 法classical :经典⽅法, 包括向前欧拉法, 改进欧拉法, 2、3、4阶龙格库塔法,Sdams-Bashford ⽅法等gear :吉尔单步法mgear :吉尔多步法2.5.1变步长龙格库塔法下⾯⽤4~5阶Runge-Kutta-Fehlberg 法求解van der Pol ⽅程:=′==+′??′′1.0)0(,0)0(0)1(2y y y y y y > ODE:=diff(y(t),t$2)-(1-y(t)^2)*diff(y(t),t)+y(t)=0;:= ODE = ? + 2t 2()y t () ? 1()y t 2t ()y t ()y t 0> initvals:=y(0)=0,D(y)(0)=-0.1;:= initvals , = ()y 00 = ()()D y 0-.1> F:=dsolve({ODE,initvals},y(t),type=numeric);:= F proc ()... end proc rkf45_x此时, 函数返回的是⼀个函数, 可以在给定的数值点上对它求值:> F(0);,, = t 0. = ()y t 0. = t ()y t -.1 > F(1);,, = t 1. = ()y t -.144768589749425608 = ??t ()y t -.178104066128215944 可以看到, F 给出的是⼀个包括t 、y(t)、D(y)(t)在内的有序表, 它对于每⼀个时间点可以给出⼀组数值表达式. 有序表的每⼀项是⼀个等式, 可对其作图描述. >plot('rhs(F(t)[2])', t=0..15, title="solution of the Van de Pol's Equation");> plots[odeplot](F,[t,y(t)],0..15,title="solution of the Van de Pol's Equation");2.5.2吉尔法求解刚性⽅程在科学和⼯程计算中, 常常会遇到这样⼀类常微分⽅程问题, 它可以表⽰成⽅程组:00)(),,(y t y y t f y ==′, 称其为刚性⽅程, 其解的分量数量相差很⼤, 分量的变化速度也相差很⼤. 如果⽤常规⽅法求解, 为了使变量有⾜够⾼的精度, 必须取很⼩的步长, ⽽为了使慢变分量达到近似的稳态解, 则需要很长的时间, 这样⽤⼩步长⼤时间跨度的计算, 必定造成庞⼤的计算量, ⽽且会使误差不断积累. 吉尔法是专门⽤来求解刚性⽅程的⼀种数值⽅法.> ODE:=diff(u(t),t)=-2000*u(t)+999.75*v(t)+1000.25,diff(v(t),t)=u(t)-v(t);:=ODE , = ??t ()u t ? + + 2000()u t 999.75()v t 1000.25 = ??t()v t ? ()u t ()v t > initvals:=u(0)=0,v(0)=-2; := initvals , = ()u 00 = ()v 0-2> ansl:=dsolve({ODE,initvals},{u(t),v(t)},type=numeric,method=gear);:= ansl proc ()... end proc x_gear> ansl(10,0);[,, = t 10.] = ()u t .989893921726687442 = ()v t .979787842765888594 > p1:=plots[odeplot] (ansl,[t,u(t)],0..20,color=red):p2:=plots[odeplot] (ansl,[t,v(t)],0..20,color=blue):plots[display] ({p1,p2}, title="Solution of a stiff equation");2.5.3 经典数值⽅法Maple中常微分⽅程数值解法中有⼀类被称作是“经典”(classical)⽅法. 当然, 称其为经典⽅法不是因为它们常⽤或是精度⾼, ⽽是因为它们的形式简单, 经常被⽤于计算⽅法课上的教学内容. 它们是⼀些常见的固定步长⽅法, 在dsolve中⽤参数method=classical[⽅法名称], 如果不特别指出, 将默认采⽤向前欧拉法. 主要有:foreuler:向前欧拉法(默认)hunform:Heun公式法(梯形⽅法, 改进欧拉法)imply:改进多项式法rk2:⼆阶龙格库塔法rk3:三阶龙格库塔法rk4:四阶龙格库塔法adambash:Adams-Bashford⽅法(预测法)abmoulton:Adams-Bashford-Moulton⽅法(预测法)下⾯给出微分⽅程数值⽅法的参数表:参数名参数类型参数⽤途参数⽤法initial 浮点数的⼀维数组指定初值向量number 正整数指定向量个数output 'procedurelist'(默认)或'listprocedure'指定⽣成单个函数或多个函数的有序表Procedurelis:单个函数, 返回有序表Listprocedure:函数的有序表procedure ⼦程序名⽤⼦程序形式指定第⼀尖常微分⽅程组的右边部分参数1:未知函数的个数参数2:⾃变量参数3:函数向量参数4:导函数向量start 浮点数⾃变量起始值startinit 布尔量(默认FALSE) 指定数值积分是否总是从起始值开始对dverk78不适⽤value 浮点数向量(⼀维数组) 指定需要输出函数值的⾃变量数值点如果给定, 结果是⼀个22×的矩阵. 元素[1,1]是⼀个向量, 含⾃变量名和函数名称; 元素[2,1]是⼀个数值矩阵, 其中第⼀列value的输⼊相同, 其他列中是相应的函数值另外, 还有⼀些特殊的附加参数:maxfun :整数类型, ⽤于最⼤的函数值数量, 默认值50000, 为负数时表⽰⽆限制 corrections :正整数类型, 指定每步修正值数量, 在abmoulton 中使⽤, 建议值≤4 stepsize :浮点数值, 指定步长下⾯看⼀个简单的例⼦:> ODE:=diff(y(x),x)=y(x)-2*x/y(x);:=ODE = ?()y x ? ()y x 2x > initvals:=y(0)=1;:= initvals = ()y 01> sol1:=dsolve({ODE,initvals},y(x),numeric,method=classical,stepsize=0.1,start=0);:= sol1proc ()... end proc x_classical⽽其解析解为:> sol2:=dsolve({diff(y(x),x)=y(x)-2*x/y(x), y(0)=1}, y(x));:= sol2 = ()y x + 2x 1将两者图形同时绘制在同⼀坐标系中⽐较, 可以发现, 在经过⼀段时间后, 欧拉法的数值结果会产⽣较⼤误差.> plot({rhs(sol2),'rhs(sol1(x)[2])'},x=0..2);求解微⽅程, ⽆论使⽤什么⽅法或者加⼊什么选项, 求解完成后必须利⽤相关数学知识进⾏逻辑判断, 绝对不对简单迷信Maple 给出的结果, 否则很有可能得到⼀个对于⽅程本⾝也许还看得过去, 但在数学或者物理意义上不合理的解.2.6摄动法求解常微分⽅程由于微分⽅程求解的复杂性, ⼀般微分⽅程常常不能求得精确解析解, 需要借助其它⽅法求得近似解或数值解, 或者两种⽅法兼⽽有之. 摄动法是重要的近似求解⽅法.摄动法⼜称⼩参数法, 它处理含⼩参数ε的系统, ⼀般当ε=0时可求得解x 0. 于是可把原系统的解展成ε的幂级数, 若这个级数当L +++=2210εεx x x x ε→0时⼀致收敛,则称正则摄动, 否则称奇异摄动. 摄动法的种类繁多, 最有代表性的是庞加莱—林斯泰特(Poicare-Lindstedt )法, 在此,我们以该⽅法求解van der Pol ⽅程:0)1(2=+′??′′y y y y ε当ε=0时该⽅程退化为数学单摆的常微分⽅程, 当ε=1时为3.5讨论的情况, 对任意ε, 该微分⽅程拥有⼀个渐进稳定的周期解, 称为极限环.由于van der Pol ⽅程中没有显式的时间依赖项, 不失⼀般性, 设初值为y(0)=0. 在庞加莱—林斯泰特法中, 时间通过变换拉伸:t ωτ=, 其中∑∞==0i i i εωω对于)(τy , van der Pol ⽅程变为:0)1(22=+′??′′y y y y ωεωrestart:diff(y(t),t$2)-epsilon*(1-y(t)^2)*diff(y(t),t)+y(t)=0;= ? +2t 2()y t ε() ? 1()y t 2()y t ()y t 0 > ODE:=DEtools[Dchangevar]({t=tau/omega,y(t)=y(tau)},%,t,tau);:= ODE = ? + ω22τ2()y τε() ? 1()y τ2ωτ()y τ()y τ0 > e_order:=6:> macro(e=epsilon,t=tau):> alias(seq(y[i]=eta[i](tau),i=0..e_order)):> e:=()->e:> for i from 0 to e_order doeta[i]:=t->eta[i](t)od:> omega:=1+sum('w[i]*e^i','i'=1..e_order);:= ω + + + + + + 1w 1εw 2ε2w 3ε3w 4ε4w 5ε5w 6ε6> y:=sum('eta[i]*e^i','i'=0..e_order);:= y + + + + + + η0η1εη2ε2η3ε3η4ε4η5ε5η6ε6> deqn:=simplify(collect(ODE,e),{e^(e_order+1)=0}):> for i from 0 to e_order doode[i]:=coeff(lhs(deqn),e,i)=0od:> ode[0];= + y 02τ2y 00> ode[1];= ? + + + 2τ2y 1τy 0y 12w 12τ2y 0??τy 0y 020 > ode[2];= + ? + + + ? + + + τy 0w 1y 022τy 0y 0y 1τy 12τy 2y 2τy 1y 02τy 0w 12w 12τy 122τy 0w 22τy 0w 120 > dsolve({ode[0],eta[0] (0)=0,D(eta[0])(0)=C[1]},eta[0](t));= y 0C 1()sin τ> eta[0]:=unapply(rhs(%),t);:= η0 →τC 1()sin τ> ode[1];= ? + ? +2τ2y 1C 1()cos τy 12w 1C 1()sin τC 13()cos τ()sin τ20 > map(combine,ode[1],'trig');= ? + ? + ? ??????2τ2y 1C 1()cos τy 12w 1C 1()sin τ14C 13()cos τ14C 13()cos 3τ0 > ode[1]:=map(collect,%,[sin(t),cos(t)]);:= ode 1 = ? + + + ? 2w 1C 1()sin τ + C 11C 13()cos τ2τ2y 1y 11C 13()cos 3τ0 > dsolve({ode[1],eta[1](0)=0,D(eta[1])(0)=C[2]},eta[1](t),method=laplace);= y 1 + ? ????????? + + 18C 1() ? C 12() + C 12τw 1C 1C 2()sin τ? 132C 13C 1τw 1()cos τ132C 13()cos 3τ > map(collect,%,[sin(t),cos(t),t]);= y 1 + ? ????????? + + 18C 1() ? C 12() + C 12τw 1C 1C 2()sin τ? 132C 13C 1τw 1()cos τ132C 13()cos 3τ> solve({coeff(lhs(ode[1]),sin(t))=0,coeff(lhs(ode[1]),cos(t))=0});,,{}, = C 10 = w 1w 1{}, = C 12 = w 10{}, = C 1-2 = w 10> w[1]:=0:C[1]:=-2:> ode[1];= + +2τ2y 1y 12()cos 3τ0 > dsolve({ode[1],eta[1](0)=0,D(eta[1])(0)=C[2]},eta[1](t),method=laplace);= y 1 ? + 14()cos 3τ14()cos τC 2()sin τ > eta[1]:=unapply(rhs(%),tau);:= η1 →τ ? + 14()cos 3τ14()cos τC 2()sin τ > map(combine,ode[2],'trig'):> ode[2]:=map(collect,%,[sin(t),sin(3*t),cos(t),cos(3*t)]);:= ode 2 = + + ? + ? + + 144w 2()sin τ54()sin 5τ2τ2y 232()sin 3τ2C 2()cos τ3C 2()cos 3τy 20 > solve({coeff(lhs(ode[2]),sin(t))=0,coeff(lhs(ode[2]),cos(t))=0}); {}, = C 20 = w 2-116> assign(%):> dsolve({ode[2],eta[2](0)=0,D(eta[2])(0)=C[3]},eta[2](t),method=laplace):> eta[2]:=unapply(rhs(%),t);:= η2 →τ? + + 316()sin 3τ+ 2996C 3()sin τ596()sin 5τ > for i from 0 to e_order domap(combine,ode[i],'trig'):ode[i]:=map(collect,%,[seq(sin((2*j+1)*t),j=0..i),seq(cos((2*j+1)*t),j=0..i)]):solve({coeff(lhs(ode[i]),sin(t))=0,coeff(lhs(ode[i]),cos(t))=0}):assign(%):dsolve({ode[i],eta[i](0)=0,D(eta[i])(0)=C[i+1]},eta[i](t),method=laplace):collect(%,[seq(sin((2*j+1)*t),j=0..i),seq(cos((2*j+1)*t),j=0..i)]):eta[i]:=unapply(rhs(%),t);od:> omega;+ + 11ε217ε435ε6。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

第四章IBM ILOG CPLEX在高速铁路列车运行图编制中的应用 数据文件:
例1:假定2列中速列车,1列高速列车,5个车站。
数据文件编码为:
nctrain=2; nhtrain=1; nstation=5; a=[4,4,4,4,4]; d=[3,3,3,3,3]; r=[[24,8,20,38],[24,8,20,38],[16,6,13,25]]; b=[[1,1,1,1,1],[1,1,1,1,1],[1,1,1,1,1]]; c=[[2,2,2,2,2],[2,2,2,2,2],[2,2,2,2,2]]; w=[[0,0,0,0,0],[0,0,0,0,0],[0,0,0,0,0]]; M=100000000; e=[[0,0,0,0,0],[0,0,0,0,0],[0,0,0,0,0]];
停站时间约束
yik xi变
k 2Sik ik qij (i 1, 2, j
, N ; k 2,3,
, m 1)
, N ; k 2,3,
, m 1)
第四章IBM ILOG CPLEX在高速铁路列车运行图编制中的应用
第四章IBM ILOG CPLEX在高速铁路列车运行图编制中的应用 例2求解结果:
// solution (optimal) with objective 1606 x = [[0 17 39 59 101 129 152] y = [[0 17 40 59 102 129 152] [0 21 43 63 106 135 159] [3 22 43 65 106 136 159] [0 26 54 79 131 164 191] [6 26 55 81 132 164 191] [0 59 81 101 143 171 194] [42 59 82 101 144 171 194] [0 63 85 105 147 174 197] [45 64 85 107 147 174 197] [0 66 88 108 151 180 204] [48 67 88 110 151 181 204] [0 73 102 127 180 215 243] [51 74 103 129 181 216 243] [0 110 131 151 194 223 246] [94 110 131 153 194 223 246] [0 114 136 156 198 226 249] [97 114 137 156 199 226 249] [0 117 140 162 206 234 257]] [100 117 141 164 207 234 257]]
第四章IBM ILOG CPLEX在高速铁路列车运行图编制中的应用
(二)最基本高速铁路运行图问题的数学模型:
目标函数:
min z [xim yi1 ]
i
N
车站间隔时间与越行约束:
l xik ak y jk M (1 qij )(k 2,3, r k
l y jk d k yik M ( qij )(k 2,3, l k 1 m
第四章IBM ILOG CPLEX在高速铁路列车运行图编制中的应用 约束: 高速铁路列车运行图根据区段和线路实际情况的不同,需要 考虑如下约束: 1. 车站间隔时间与越行约束。 2. 区间运行时间约束。 3. 停站时间约束。 4. 天窗时间约束。 5. 到发线约束 6. 机车数量约束 7.其他约束 其中,约束1、2、3称为最基本的约束,此外,根据实际情况 添加相应的约束,相应的模型的复杂程度会有所增加。
第四章IBM ILOG CPLEX在高速铁路列车运行图编制中的应用
(一)高速铁路运行图问题的描述
高铁列车运行图指在铁路网上同时考虑各客运车站衔接区 间各列车安全或作业时间间隔、区间运行时间、起停附加 时间、天窗起止时间和最少停站时间的作业时间等标准的 约束,确定列车开行方案中各次列车在各个车站出发、到 达或通过时间,还要确保列车在合理时间域内始发和终到, 在这些约束下使得各次列车总旅行时间最短,并且保证运 行线间的鲁棒性,得到方案中各次列车运行轨迹的时间— 空间图解。
(三) 计算机模型的OPL语言:
模型文件:
int nctrain=...;//200km/h动车组数量 int nhtrain=...;//300 km/h动车组数量 int nstation=...;//车站数 range station=1..nstation; range section=1..nstation-1; //range ctrain=1..nctrain; range htrain=1..nctrain+nhtrain; int a[station]=...; //车站高速列车到达间隔 int d[station]=...;//车站高速列车发车间隔
第四章IBM ILOG CPLEX在高速铁路列车运行图编制中的应用 例1.某高速列车运行图,如下图所示,有7列快车3列车慢车:
714.72
602.4
468.98
269.25
180.69
78.45
10
20
30
40
50
60
70
80
90
100 110 120 130 140 150 160 170 180 190 200 210 220 230 240
第四章IBM ILOG CPLEX在高速铁路列车运行图编制中的应用 模型文件(续):
int r[htrain][section]=...;//不同列车在相应区间运行时间 int b[htrain][station]=...;//列车在车站停车附件时分 int c[htrain][station]=...;//列车在车站起车附件时分 int w[htrain][station]=...;//列车在车站的停站时间 int M=...; //int e[htrain][station]=...; dvar int+ x[htrain][station] in 0..1000;//列车到达车站时分,整数变量 dvar int+ y[htrain][station]in 0..1000;//列车车站出发时分,整数变量 dvar int+ q[htrain][htrain][station] in 0..1;//0-1变量,两列车是否在车站越 行,是取1,否为0 dvar int+ s[htrain][station] in 0..1;// 0-1变量,列车是否在车站要求停站, 是取1.否为0.
第四章IBM ILOG CPLEX在高速铁路列车运行图编制中的应用 例1求解结果:
// solution (optimal) with objective 240 x = [[0 24 32 52 90] [0 28 36 56 94] [0 54 60 73 98]] y = [[0 24 32 52 90]
[4 28 36 56 94]
[38 54 60 73 98]]
第四章IBM ILOG CPLEX在高速铁路列车运行图编制中的应用
例2:假定3列中速列车,7列高速列车,7个车站。 数据文件编码为:
/*中速列车数量*/ nctrain=3; /*高速列车数量*/ nhtrain=7; /*车站数量*/ nstation=7; /*发车间隔时间*/ a=[3,3,3,3,3,3,3]; /*到达间隔时间*/ d=[3,3,3,3,3,3,3]; /*区间运行时间*/ r=[[16,20,18,40,27,23],[16,20,18,40,27,23],[19,25,21,48,32,27],[16,20,18, 40,27,23],[16,20,18,40,27,23],[16,20,18,40,27,23],[19,25,21,48,32,27],[1 6,20,18,40,27,23],[16,20,18,40,27,23],[16,20,18,40,27,23]];
, m 1)
, m 1)
i, j 1, 2,
N,i j
第四章IBM ILOG CPLEX在高速铁路列车运行图编制中的应用
区间运行时间约束:
xik 1 xik rik ik Sik ik 1Sik 1 (i 1, 2, , N ; k 1, 2, , m 1)
第四章IBM ILOG CPLEX在高速铁路列车运行图编制中的应用
数据文件编码(续): /*停站附加时间*/ b=[[1,1,1,1,1,1,1],[1,1,1,1,1,1,1],[1,1,1,1,1,1,1],[1,1,1,1,1,1,1],[1,1,1,1,1,1, 1],[1,1,1,1,1,1,1],[1,1,1,1,1,1,1],[1,1,1,1,1,1,1],[1,1,1,1,1,1,1],[1,1,1,1,1,1, 1]]; /*启动附加时间*/ c=[[2,2,2,2,2,2,2],[2,2,2,2,2,2,2],[2,2,2,2,2,2,2],[2,2,2,2,2,2,2],[2,2,2,2,2,2, 2],[2,2,2,2,2,2,2],[2,2,2,2,2,2,2],[2,2,2,2,2,2,2],[2,2,2,2,2,2,2],[2,2,2,2,2,2, 2]]; /*停站时间矩阵*/ w=[[0,0,1,0,1,0,0],[0,1,0,2,0,1,0],[0,0,1,2,1,0,0],[0,0,1,0,1,0,0],[0,0,0,2,0,0, 0],[0,1,0,2,0,1,0],[0,1,1,2,1,1,0],[0,0,0,2,0,0,0],[0,0,1,0,1,0,0],[0,0,1,2,1,0, 0]]; M=100000000;
相关文档
最新文档