当前位置： >>
Matlab与化工数值计算
精通MATLAB语言，有MATLAB编程问题的朋友，请直接联系我！ 联系方式： QQ： Email： MATLAB博客：http://top99./_d.html如果我QQ不在线，可以将问题发到我的邮箱，或者在我的MATLAB 博客留言，第一时间答复你！简单问题请
直接留言，复杂问题可以 提供解决思路！同探讨，共进步！还可以为公司和科研单位设计各类算法，欢迎前来咨询！郑重声明： 本资源来源于网络， 仅限交流使用， 请勿用于商业用途！Matlab与化工数值计算第一讲 简介与基本数学运算隋志军 化工学院软件应用教科组 联系方式：zhjsui@
实验16楼605室化学工程师的任务工厂运行 系统特性 过程设计 数学模型 反应特征 过程特性 设备特性 工艺开发 数学模型化学工程专业数学模型类型非线性方程（组）1 f = 4.0 0.4 lg(Re gen f 1? n ' 2 ) ? 1.2 n' 0.75 n'常微分方程（组）f & '+ ff &+ β (1 ? ( f ' ) 2 ) = 0偏微分方程（组）?u1 ? 2 u1 = 0.024 2 ? F (u1 ? u 2 ) ?t ?x非线性模型，难以获得 解析解，必须采用数值 解法 模型的数值解法是应用 数学的一个分支，通常 称为计算数学（数值分 析，数值方法）? 2u 2 ?u 2 + F (u 1 ? u 2 ) = 0.170 2 ?t ?x化学工程常用软件数学软件：? Matlab ? Mathematica ? Mathcad ? Maple ? Staticstica化工模拟软件：PRO/II (SimSci) AspenPlus ChemCAD Flowtran Superpro Designer Fluent CHEMKIN本课程的学习目的学会Matlab的使用，可以利 用Matlab求解较为复杂的化 工数学模型 Matlab 对于数值分析 的内容不过多 涉及，只注意 数值计算结果 的准确性 化工专业知识 作为背景，不 涉及模型的推 导，注重模型 求解过程的方 法与技巧数值分析化工数学模型本课程基本内容 ? ? ? ? ? ? ? 第一讲 Matlab简介与基本数学运算 第二讲 非线性方程组求解与迭代法 第三讲 矩阵操作与线性方程组求解 第四讲 插值、拟合与数值微分、积分 第五讲 常微分方程数值解 第六讲 偏微分方程数值解 第七讲 统计初步与最优化方法学习本课程的注意事项? 学好本课程的唯一途径是多上机实践 ? 数值计算效率和效果的保证有很多技巧,可以参考数 值方法（数值分析）方面的教科书 刘则毅，科学计算技术与Matlab，科学出版社 同济大学计算数学教研室，现代数值数学和计 算，同济大学出版社 黄华江，实用化工计算机模拟，化学工业出版社 张志涌，精通Matlab6.5版，北京航空航天大学出 版社 ? 对于数值计算的结果，应注意分析结果的意义Matlab简介Matlab是Matrix Labotary的 缩写，最初是美国新墨西哥 大学Moler教授编写的 LINPACK和EISPACK接口 程序1984年，MathWorks公司创 建，MATLAB正式推向市场 20世纪90年代以来，MATLAB已 成为数值计算软件的佼佼者Prof. Cleve Moler Jack LittleMatlab简介? MATLAB具有用法简单、 灵活、结构性强、延展性 好等优点，逐渐成为科技 计算、视图交互系统和程 序中的首选语言工具。C 功能强大的数值运算功 能 C 强大的图形处理能力 C 高级但简单的程序环境 C 丰富的工具箱与模块集 C 易于扩充开始的问题计算在1/2英寸不锈钢管中，以2000lb/hr流量输送水，当水的 温度为10、20、30、40、50、60、70、80℃时，压降分别为 多少？ M μ ΔP = 牛顿流体在不锈钢管中的流动压降可由下式估算： 20000D ρ1. 8 0. 2 4. 8其中，摩擦压降，psi/(100英尺等量管长)；M，质量流量，lb/hr；μ，粘 度，cP；ρ，密度，lb/ft3，D，管径，inch。流体密度可由下式描述： = A ? B ? (1?T T ρC)nρ，g/ml；对于水，A＝0.34710；B＝0.2740；Tc＝647.13K；n＝0.28571。μ 流体粘度由下式描述：log 10= A + B / T + CT + DT 2μ，cP；对于水，A=-10.2158；B=1.7925E3；C＝1.7730E-2；D=-1.2631E-05。Matlab窗口介绍变量空间 命令窗口 当前路径命令历史Matlab的通用命令命令 cd type clf clc hold path load quit whos 命令说明 显示或改变工作目录 显示文件内容 清除图形窗口 清除命令窗口内容 图形保持开关 显示搜索目录 加载指定文件变量 退出Matlab 变量查看 命令 dir clear pack echo disp save diary ！ 命令说明 显示目录文件 清除内存变量 收集内存碎片，扩大 内存空间 命令窗口信息显示开 关 显示变量或文字内容 保存内存变量到指定 文件 日志文件命令 调用DOS命令通过Help学习Matlab在命令窗口中键入&& help，则显示以下内容：matlab\general matlab\ops matlab\lang matlab\elmat matlab\elfun matlab\specfun matlab\matfun matlab\datafun matlab\polyfun matlab\funfun matlab\sparfun matlab\scribe matlab\graph2d matlab\graph3d ……… - General purpose commands. - Operators and special characters. - Programming language constructs. - Elementary matrices and matrix manipulation. - Elementary math functions. - Specialized math functions. - Matrix functions - numerical linear algebra. - Data analysis and Fourier transforms. - Interpolation and polynomials. - Function functions and ODE solvers. - Sparse matrices. - Annotation and Plot Editing. - Two dimensional graphs. - Three dimensional graphs.Help＋主题名称&& help ops Operators and special characters. Arithmetic operators. plus - Plus uplus - Unary plus minus - Minus uminus - Unary minus mtimes - Matrix multiply times - Array multiply mpower - Matrix power power - Array power mldivide - Backslash or left matrix divide mrdivide - Slash or right matrix divide ldivide - Left array divide rdivide - Right array divide+ + * .* ^ .^ \ / .\ ./基本算术运算符运 加 矩阵乘 矩阵左除 矩阵右除 幂次方 算 符 + * \ / ^ 号 运 算 减 数组相乘 数组左除 数组右除 数组幂次方 符 .* .\ ./ .^ 号Help＋函数名&& help power .^ Array power. Z = X.^Y denotes element-by-element powers. X and Y must have the same dimensions unless one is a scalar. A scalar can operate into anything. C = POWER(A,B) is called for the syntax 'A .^ B' when A or B is an object. Help+函数名可获得详细的函数使用方法Matlab语言的标点标点 ： ； ， （） [] {} 定义 向量和矩阵的多 种功能 区分行及取消行 显示 区分列及函数参 数分隔符 指定运算过程的 次序等 矩阵定义 构成单元数组 标点 . ... ％ ！ ＝ ‘ 定义 小数点及结构体 域的访问 续行符 注释符，百分号 调用dos操作命 令 赋值标记 字符串标示符数值类型?分类方法一C双精度型 （系统默认类型） C单精度型 C带符号整数 C无符号整数?分类方法二C标量 C向量 C数组?分类方法三C实数 C复数数值的表示 以下表达方式均合法： 345 -99 1.3e-3 4.5e33 [1 2 3] [1;2;3] [1 2; 2 11] 3+3i 6-8j0.01基本数学运算符的使用计算以下表达式的值： 1) [1 2 3]*[3 2 1] 2) [1 2 3].*[3 2 1] 3) [1 2 3]^2 4) [1 2 3].^2 5) 1+3*2^2 6) (3*2)^2 7) (3*2)^2; 8) (-8)^(1/3)Matlab的计算器功能&& *(10^(-10.5e3/283+1.773e-2*e- 5*283^2))^0.2/(^4.8*(0.^(-()^0.28574))/0.2323) 回车可以得到结果 ans = 287.8245命令的窗口的快捷键快捷键 ↑，Crtl+P ↓，Crtl+N ←，Crtl+B →，Crtl+F Crtl+← Ctrl+C 作用 回调上一行 回调下一行 回移上一字符 前移下一字符 左移一单词 终止正在运行的 程序 快捷键 Crtl+→ Crtl+A，Home Crtl+E，End Crtl+U，Esc Crtl+K 作用 右移一单词 移至行首 移至行末 删除一行 从光标删除至行 末数学函数（elfun）类型三角函数函sin(x) asin(x) cos(x)数正弦值 反正弦值 余弦值 反余弦值 正切 指数运算 自然对数 求平方根 求绝对值 取出复数的虚部 取出复数的实部 复数共轭 四舍五入 求余数 整数x和y的最小公倍数 整数x和y的最大公约数含义acos(x) tan(x) 指数函数 exp(x) log(x) sqrt(x) 复数函数 abs(x) imag(x) real(x) conj(x) 数论函数 round(x) mod(x,y) lcm(x,y) gcd(x,y)基本数学函数的使用 计算以下表达式的值： 1) sin(30) 2) sind(30) 3) exp([1 2 3]) 4) log10([10 100 1000]) 5) abs(3+4i) 6) abs(-5)format命令 MATLAB命令format short format short e format long format long e format rat format hex format bank含短格式义3.1416范3.例短格式科学格式 长格式 长格式科学格式 有理格式 十六进制格式 银行格式3.79 3.793e+000 355/113 442d18 3.14程序的组成变量 数据输入键盘输入 文件输入变量 运算数学运算 关系运算 逻辑运算 流程控制数据输出屏幕输出 文件输出 图形输出变量? 变量的命名方式：C 变量名由字母、数字和下划线组成； C 变量名中的英文字母大小写是有区别的； C 变量名的最大长度是有规定的? 不同版本的系统规定不同：19个字符、31或63个 字符等 ? 可调用namelengthmax函数得到系统规定长度变量的使用&&clear %删除工作区中所有定义过的变量 &&whos %查看当前工作区内变量信息，无显示表示没 有定义的变量 && xy=1; yx=2; %对变量赋值 && xy %查看变量xy的当前数值 xy = 1 && whos Name Size Bytes Class xy 1x1 8 double array yx 1x1 8 double array Grand total is 2 elements using 16 bytes && clear xy yx %删除变量xy及yx && whos && xy %这时变量xy已经不存在了 ??? Undefined function or variable 'xy'.MATLAB系统的特殊变量和常数特殊变量ans pi inf或Inf eps realmax realmin NaN或nan i或j nargin nargout lasterr lastwarn意圆周率π（= 3.1415926...） 无穷大∞值 浮点运算的相对精度2^(-52) 最大的正浮点数，2^(1024)-1 最小的正浮点数，2^(-1022) 不定量 虚数单位 函数输入参数个数 函数输出参数个数 存放最新的错误信息 存放最新的警告信息义如果用户未定义变量名，系统用于计算结果存储的默认变量名MATLAB数据类型? ? ? ? ? 数值（标量，向量，数组） 字符串 单元数组（cell array） 结构体（structure） 函数句柄向量的生成1) 直接输入向量格式上要求向量元素需要用“[]”括起来，元素之间可以用空格、逗号或分号分隔。 用空格和逗号分隔生成行向量，用分号生成列向量。2) 利用冒号生成向量冒号表达式的基本形式为：x＝x0:step:xn 若step＝1，则此项输入可以忽略。3) linspace函数可以使用linspace函数生成线性等分向量： y＝linspace(x1,x2) 生成（1*100）维行向量，y(1)=x1，y(100)=x2 y＝linspace(x1,x2, n) 生成（1*n）维行向量，y(1)=x1，y(n)=x24) logspace函数logspace用于生成对数等分向量，格式如下： y＝logspace(x1,x2,n) 生成（1*n）维对数等分向量，y(1)=10^x1， y(n)=10^x2；n可以省略，此时其默认值为50。向量的运算1) 向量加减与数加减向量的加减与数加减的形式与普通标量加减相同2) 向量的点积、叉积与混合积的实现点积：向量的点积由函数dot实现。dot(a,b)返回向量a和b的数量点积，其 中a，b必须同维。 叉积：叉积由cross函数实现。向量a,b必须为三维向量 混合积：可由以下命令实现，dot(a,cross(b,c))3)向量的数乘、数组乘和向量乘例：当a＝[1:1:3]; b=[2:2:6]时，以下命令的运行结果是什么？ 1)&&a1=2*a 2) && a2=a.*b 3) &&a3=a*b字符串类型 ? 字符串：包含在一对单引号中的字符集合&& s='hello, MATLAB's= hello, MATLAB && whos Name Size Bytes Class s 1x13 26 char array Grand total is 13 elements using 26 bytes %定义字符串变量s单元数组（Cell Array）单元数组是MATLAB数组的一种特殊数据类型，它用于保存 不同类型和/或不同大小的数据。 三种直接赋值方式1. 单元下标用括号“（）”括起来，而单元的内容用“{}”括起来，如： &&clear all &&a(1,1)={[1 2;3 4]}; &&a(1,2)={[0 1]}; &&a(2,1)={'Hello'}; &&a(2,2)={2+3i}2.单元下标用“{}”括起来，而赋值语句等式右边的单元内容用“[]”括起来：&&a{1,1}=[1 2;3 4]; &&a{1,2}=[0 1]; &&a{2,1}=’Hello’; %右边只有一个元素时可省略去“[]” &&a{2,2}=2+3i3. 直接使用{}&&a={[1 2;3 4],[0 1],’Hello’,2+3i}单元数组的操作显示单元数组的命令&&a % 显示单元数组a的信息 % 显示单元数组a的完整内容&&celldisp(a)先使用函数cell()创建空的单元数组，然后再赋值：&&b=cell(2,3) 赋值方法同直接赋值方式。对单元数组元素的操作&&c=a{1,2} “()”。 % 将单元数组a的{1,2}元素赋给变量c，注意是“{}”,而不是结构体? 与C语言类似，MATLAB结构体用于存取相关的 数据，它由一组称为域（fields）的成员变量 （向量）构成，每一个域可以为不同的MATLAB 数据类型。 ? 结构数组的定义有两种方法，一种是直接赋值， 另一种是使用strct()函数。结构体的赋值&&student.name=’Zhang Jun’; &&student.major=’Chemical Engineering’; &&student.subject=[‘英语 ’,’政治 ’,’数学 ’,’化工原理 ’,’物理化学 ’]; &&student.entrance_exam=[62 68 72 82 90];&&student(2).name=’Li Xia’; &&student(2).major=’Chemical Engineering’; &&student(2).subject=[‘英语 ’,’政治 ’,’数学 ’,’化工原理 ’,’物理化学 ’]; &&student(2).entrance_exam=[60 72 68 85 88];struct_array_name=structure(‘field1’,values1,’field2’,values2,…) Student=struct(‘name’,’Zhang Jun’, ‘major’,’Chemical Engineering’)管道压降的计算T=[283:10:353]; M=2000;D=0.5; density.A=0.3471;density.B=0.274;density.Tc=647.13;density.n=0.28571; Rho=(density.A.*density.B.^(-(1-T./density.Tc).^density.n))/0.2323; mu.A=-10.2158;mu.B=1.7925e3;mu.C=1.773e-2;mu.D=-1.2631e-5; mu=10.^(mu.A+mu.B./T+mu.C.*T+mu.D.*T.^2); deltP=(M^1.8)*(mu.^0.2)./(20000*D^4.8.*Rho)函数文件和Script文件Script文件Script仅仅是一连串可执行的MATLAB命令，它具有全局性 Script文件中不能定义函数函数函数定义的一般格式： function [y1,y2,…,yn] = FuncName(x1,x2,…,xn) %函数声明语句 y1 = … % （表达式1） y2 = … % （表达式2） … yn = … % （表达式n） 其中，输入参数为x1,x2,…,xn,输出参数为y1,y2,…,yn。各参数可以是 标量、向量或矩阵。脚本编辑窗口函数文件的编写编写一个函数，计算本章开始问题中流体的粘度， 函数要求输出粘度的计算值： 1. 函数声明语句：function vis=viscosity() 2. 变量的传递 1. 通过调用函数传递 ? vis=viscosity(A,B,C,D,T) 2. 通过全局变量传递 ? 利用global命令，在主函数和子函数中予以声明 3. 编写表达式函数的调用1. 在调用函数的主函数中，直接采用函数名调用 2. 通过函数句柄调用管道压降的计算函数function Cha1demo4_7 global mu T T=[283:10:353]; M=2000;D=0.5; density.A=0.3471; density.B=0.2740; density.Tc=647.13; density.n=0.28571; Rho=(density.A.*density.B.^(-(1-T./density.Tc).^density.n))/0.2323 mu.A=-10.2158;mu.B=1.7925e3;mu.C=1.773e-2;mu.D=-1.2631e-5; Mu=viscosity deltP=(M^1.8)*(Mu.^0.2)./(20000*D^4.8.*Rho) %---------------------------------------------------------------function vis=viscosity global mu T vis=10.^(mu.A+mu.B./T+mu.C.*T+mu.D.*T.^2);函数句柄? 创建一个函数句柄，可用于保存函数的所有信 息，以便将来对它进行调用。 ? 函数句柄可以作为参数传递给其他函数，或与 feval函数一起使用，以调用该函数句柄所属的函 数。 ? 使用函数句柄还可以减少定义函数的文件个数， 改善重复操作的性能，保证函数计算的可靠性。 ? funhandle = @function_name %function_name为用户指定的函数名管道压降的计算函数function Cha1demo4_5 T=[283:10:353]; M=2000;D=0.5; density.A=0.3471; density.B=0.274 density.Tc=647.13; density.n=0.28571; Rho=(density.A.*density.B.^(-(1-T./density.Tc).^density.n))/0.2323 mu.A=-10.2158;mu.B=1.7925e3;mu.C=1.773e-2;mu.D=-1.2631e-5; Mu=feval(@viscosity,mu,T)%Mu=viscosity(mu,T) deltP=(M^1.8)*(Mu.^0.2)./(20000*D^4.8.*Rho) %---------------------------------------------------------------function vis=viscosity(mu,T) vis=10.^(mu.A+mu.B./T+mu.C.*T+mu.D.*T.^2);内联函数（inline function）? 内联函数是Matlab提供的一个对象，它的表现和函数文件 一样，但内联函数的创建比较容易 ? 内联函数的创建 inline(‘CE’) inline(‘CE’,arg1,arg2,...) inline(‘CE’,n) ? 涉及内联函数性质的指令 class(inline_fun) 内联函数类型 char(inline_fun) 给出内联函数计算公式 argnames(inline_fun) 给出内联函数的输入变量 vectorize(inline_fun) 使内联函数适用于数组运算 ? Matlab中许多“泛函”函数都是采用inline，从而具备了适应 各种被处理函数形式的能力内联函数的应用F1=inline('sin(rho)/rho') f1=F1(2)FF1=vectorize(F1) xx=[0.5,1,1.5,2];ff1=FF1(xx)G2=inline('a*exp(x(1))*cos(x(2))','a','x') g1=G2(2,[-1,pi/3])匿名函数（anonymous function）? 匿名函数用于在命令行、函数文件或script文件中创建简 单形式的函数，避免另外定义新的函数 ? 匿名函数的定义形式 f=@(arglist)expressionf=@(x) x.^2 a=f(5) 结果：a＝25 f=@(x) x.^2; g=@(x) 3*x; h=@(x) g(f(x)); h(3) 结果：ans＝27alpha=0.9; f=@(x) sin(alpha*x)； f(pi) 结果：ans＝0.3090数据输入和输出? 数据输入C利用M文件产生数据文件 C用Load命令从MAT文件或文本文件读取数据 C用fscanf函数 C用提示输入函数input Cdlmread, importdata, xlsread函数? 数据输出C用Save命令 C用fprintf函数 C用函数disp()将结果输出至屏幕 Cdlmwirte,xlswrite函数 C 图形输出Matlab二维图形1 数据准备： 选定所要表现的范围 产生自变量采样向量 计算相应的函数值向量 选定图形窗及子图位置 缺省时，打开Figure No.1，或当前窗， 当前子图 可用指令指定图形窗号和子图号 调用（高层）绘图指令；在指令中设置 线型、色彩、数据点型 设置轴的范围与刻度、坐标分格线 图形注释： 图名、坐标名、图例、文字说明 t=pi*(0:100)/100; y=sin(t).*sin(9*t);2figure(1) ％指定1号图形窗 subplot(2,2,3) ％指定一个具有2行2列子 图图形窗中的3号子图 plot(t,y,'b-') %用蓝色实线画曲线3 4 5axis([0,pi,-1,1]) %设置轴的范围 grid on ％画坐标分格线 title('调制波形') %图名 xlabel('t');ylabel('y') %轴名 legend('sin(t)','sin(t)sin(9t)') %图例 text(2,0.5,'y=sin(t)sin(9t)') %文字说明 set(h,'MarkerSize',10) %设置数据点大小6图形的精细修饰（图柄操作）： 利用对象属性值进行设置 利用图形窗工具条进行函数Plot基本调用格式1) plot(X,'s') ? X为实向量时，以该向量元素的下标为横坐标、元素值为纵坐标画一条连续曲 线 ? X是实矩阵时，则按列绘制每列元素值相对其下标的曲线，曲线数目等于X的 列数 ? X是复数矩阵时，则按列分别以元素的实部和虚部为横、纵坐标绘制多条曲线 ? 's'是用来控制线型、色彩、数据点型的选项字符串。s可以缺省，此时曲线按 Matlab默认设置绘制。s的取值见下节 2) plot(X,Y,'s') ? X、Y是同维向量时，绘制以X、Y为横、纵坐标的曲线 ? X是向量，Y是有一维与X同维的矩阵时，则绘出多根不同色彩的曲线。曲线 数等于Y的另一维，X作为这些曲线共同的横坐标 ? X是矩阵，Y是向量时，情况与上相同，只是曲线都以Y为共同纵坐标 ? X、Y是同维矩阵时，则以X、Y对应列元素为横、纵坐标分别绘制曲线，曲线 条数等于矩阵的列数。 ? 's'的意义，与上相同。 3) plot(X1,Y1,'s1',X2,Y2,'s2',…) ? 此格式中，每个绘线“三元组”（X,Y,'s'）的结构和作用，与上相同。不同“三元 组”之间没有约束关系。曲线的色彩、线型和数据点型貌符号 线型 含义 符号 含义 实线 b 蓝色 g 绿色 虚线 r 红色 c 青色 ： 点划 线 m 品红 y 黄色 -. 双划 线 k 黑色 w 白色 --色彩符号 . + * ^ & & v含义 实心黑点 十字 八线符 上三角 右三角 左三角 下三角符号 d h o p s x含义 菱形 六角星符 空心圆圈 五角星符 方块 叉字符坐标、刻度和分格线控制指令 axis auto axis manual axis off axis on axis ij axis xy axis(V) V=[x1,x2,y1,y2 ] V=[x1,x2,y1,y2, z1,z2] 含义 缺省设置 当前坐标范围不变 取消轴背景 使用轴背景 矩阵式坐标，原点在左上 方 普通直角坐标 人工设定坐标范围。设定 值：二维，4个；三维，6 个 指令 axis equal axis fill axis image axis normal axis square axis tight axis vis3d 含义 纵横轴采用等长刻度 manual方式下起作用，使 坐标充满整个绘图区 纵、横轴采用等长刻度， 坐标框紧贴数据范围 缺省矩形坐标系 产生正方形坐标系 把数据范围设为坐标范围 保持高宽比不变图形标识图形标识命令 图名（title） 坐标轴（label） 图形注释（text） 图例（legend） 调用格式TITLE(‘text’) TITLE('text','Property1',PropertyValue1,'Property2',P ropertyValue2,...) XLABEL('text’) XLABEL('text','Property1',PropertyValue1,'Property2',Pr opertyValue2,...) TEXT(X,Y,'string')LEGEND(string1,string2,string3, ...)例题T=[283:10:353]; M=2000;D=0.5; density.A=0.3471;density.B=0.274;density.Tc=647.13;density.n=0.28571; Rho=(density.A.*density.B.^(-(1-T./density.Tc).^density.n))/0.2323; mu.A=-10.2158;mu.B=1.7925e3;mu.C=1.773e-2;mu.D=-1.2631e-5; mu=10.^(mu.A+mu.B./T+mu.C.*T+mu.D.*T.^2); deltP=(M^1.8)*(mu.^0.2)./(20000*D^4.8.*Rho)plot(T,deltP,'b-o') title('The pressure drop vs Temperature') xlabel('Temperature(^oC)') ylabel('Pressure drop (psi/equivalent feet of pipe)')图形标识的精细控制句柄图形的结构Root FigureUicontrolAxesUimenuLineSurfaceRectanglePatchImageTextLight句柄图形的访问与操作函数名 gca gcbf gcbo gcf gco copyobj delete findobj get set 说明 获得当前坐标轴对象的句柄 获得当前正在执行调用的图形对象句柄 获得当前正在执行调用的对象的句柄 获得当前图形对象的句柄 获得当前对象的句柄 复制图形对象及相应子对象 删除图形对象 依属性值查找图形对象 获得对象属性 设置对象属性例题 绘制二阶系统阶跃响应曲线例题t=6*pi*(0:100)/100;y=1-exp(-0.3*t).*cos(0.7*t); tt=t(find(abs(y-1)&0.05)); %寻找大于0.05的元素 ts=max(tt); ％寻找tt中最大的元素 plot(t,y,'r-','LineWidth',3) axis([-inf,6*pi,0.6,inf]) set(gca,'Xtick',[2*pi,4*pi,6*pi],'Ytick',[0.95,1,1.05,max(y)]) grid on title('\it y = 1 - e^{ -\alphat}cos{\omegat}') text(13.5,1.2,'\fontsize{12}{\alpha}=0.3') text(13.5,1.1,'\fontsize{12}{\omega}=0.7')plot(ts,0.95,'bo','MarkerSize',10);hold off cell_string{1}='\fontsize{12}\uparrow'; cell_string{2}='\fontsize{16} \fontname{隶书}镇定时间'; cell_string{3}='\fontsize{6} '; cell_string{4}=['\fontsize{14}\rmt_{s} = ' num2str(ts)]; text(ts,0.85,cell_string) xlabel('\fontsize{14} \bft \rightarrow') ylabel('\fontsize{14} \bfy \rightarrow')多次叠绘、双纵坐标和多子图 1)多次叠绘：hold on，使当前轴及图形保持而不被刷新；hold off，不保持当前轴及图形。2) 双坐标图：plotyy(X1,Y1,X2,Y2) 以左右不同纵轴分别绘制X1-Y1、X2-Y2两条曲线 plotyy(X1,Y1,X2,Y2,Fun) 以左右不同纵轴把X1-Y1、X2-Y2绘制成Fun指定的形式 的两条曲线3)多子图：采用subplot(m,n,k)使（m×n）幅子图中的第k个成为当前子图，再采用其它的图形 绘制指令则可将图形绘制到指定的子图中。子图序号的编制原则是：左上方为第1 幅，向右向下依次增大。双坐标曲线绘制方法画出函数 y = x sin x 曲线 和积分 s =5∫x0( x sin x ) dx在区间 [0,4] 上的4s= ∫ xxsinxdx00y=xsinx2-500.511.522.533.50 4 dx=0.1;x=0:dx:4;y=x.*sin(x);s=cumtrapz(y)* plotyy(x,y,x,s),text(0.5,0,'\fontsize{14}\ity=xsinx') sint='{\fontsize{16}\int_{\fontsize{8}0}^{ x}}'; text(2.5,3.5,['\fontsize{14}\its=',sint,'\fontsize{14}\itxsinxdx']) %梯形法求累计积分Matlab三维图形1) 三维曲线绘制命令 plot3 2) 三维网格图形绘制命令 mesh 3) 三维曲面绘制名利 surf[X0 Y0 Z0]=sphere(30); X=2*X0;Y=2*Y0;Z=2*Z0; surf(X0,Y0,Z0); shading interp hold on mesh(X,Y,Z) colormap(hot) hidden off hold off axis equal axis offMatlab图形绘制函数Matlab图形绘制函数分属于以下帮助主题 ? graph2d ? graph3d ? specgraph程序的组成Script文件或函数文件变量 数据输入键盘输入 文件输入变量 运算数学运算 关系运算 逻辑运算 流程控制数据输出屏幕输出 文件输出 图形输出MATLAB关系运算符关系运算符： 运 算 符 & == &= 号 运 算 符 & ~= &= 号大于 等于 大于等于小于 不等于 小于等于MATLAB逻辑、关系运算的规定? 逻辑和关系运算均可以在任意具有相同维数的数组之 间进行 ? 标量可以和任意维数的数组运算 ? 在所有关系表达式和逻辑表达式中，作为输入的任何 非0数均被看作是“逻辑真”，只有0是“逻辑假” ? 所有关系表达式和逻辑表达式的运算结果是一个由0 和1组成的逻辑数组 ? 逻辑数组是“数值数组”的子类，可以按数值数组实施 操作；它又有不同与普通数值数组的特性，即表示着 对对象判断的真假；可用于数组寻访等特殊场合MATLAB逻辑运算符逻辑运算符： 运 与 非 A=[0 1 1 0 ]; 算 符 & ~ 号 运 或 异或 算 符 | xor 号B=[1 1 0 0 ]A&B=[0 1 0 0] A|B=[1 1 1 0] ～A=[1 0 0 1] xor(A,B)=[1 0 1 0] 先决逻辑运算符： 先决与 && 先决或 ||逻辑、关系运算示例A=-3:3; B=3:-1:-3; L1=~(A&0) L2=~A&0 L3=~A L4=A&-2&A&1 L5=A==B&L1运算优先级：1) 2) 3) 4) 5) 6) 7) 8) 括号 逻辑否 乘除 加减 关系运算 逻辑运算 先决与 先决否高低运行结果：A = -3 -2 -1 0 1 2 3 B = 3 2 1 0 -1 -2 -3 L1 = 1 1 1 1 0 0 0 L2 = 0 0 0 1 0 0 0 L3 = 0 0 0 1 0 0 0 L4 = 0 0 1 1 0 0 0 L5 = 0 0 0 1 0 0 0MATLAB关系运算函数isempty isequal any all find isscalar isvector isnan isinf isfinite 数组是否为空 两个数组是否相等 数组有非零元素则结果为1 数组元素全非零则结果为1 数组非零元素的下标 是否为标量 是否为向量 是否为非数 是否为无穷 是否为有限for循环结构for 循环变量 = 表达式1（初值）：表达式2（步长）：表达式3（终值） statements （语句组） end字符串：包含在一对单引号中的字符集合 for i=1:10 x(i)=i; end x ? ? 为了得到高效代码，应尽量提高代码的向量化程度，避免使用循环 结构 为了得到高效代码，在循环指令之前应尽量对数组进行预定义while循环结构while cndition（表达式） statements（执行语句组） end Fibonacci数组的元素满足Fibonacci规则： ak+2=ak+ak+1，（k=1,2,...);且a1=a2=1。求该数组中第一 个大于10000的元素。a(1)=1;a(2)=1;i=2; while a(i)&=10000 a(i+1)=a(i-1)+a(i); i=i+1; end i, a(i)if-else-end分支结构if语句的一般格式： if condition1 statements %如果condition1的值为True，则执行该语句组 elseif condition2 statements %如果condition2的值为True，则执行该语句组 else statements %如果condition1和condition2的都为False，则执行该语句组 end 用for循环指令寻求Fibonacci数组中第一个大于10000的元素 n=100;a=1:100; 可以用break语句强制终止循环的运 for i=3:n 行 a(i)=a(i-1)+a(i-2); if a(i)&=10000 a(i), end end iswitch-case结构switch-case的一般格式： switch test_expr %测试表达式test_expr可以是标量或字符串 case value statements ％当test_expr值是value时，执行该语句组 case {value1,value2,…} statements ％当test_expr值是value1或value2或……时，执行该语句组 otherwise, statements endswitch…case…otherwise语句的能力与if…else…end语句类似， 但对多重选择的情况switch语句使代码更加易读。try-catch结构try-catch的一般格式： trystatements ％此组语句总被执行。若正确则跳出此结构 ％当上组语句出现执行错误后，该组语句被执行catchstatementsend?当两组语句都出错后，Matlab将 跳出该结构 ?可以采用lasterr函数查询出错原因N=4;A=[1 2 3]; try A_N=A(N) catch A_N=A(end) end lasterrA_N = 3 ans = Index exceeds matrix dimensions.本讲小结1)Matlab的基本数学运算符和运算函数的使用? 注意区别矩阵和数组的乘、除、乘方运算2)Matlab数据输入输出功能，尤其是绘图功能的实现 3)Matlab函数文件的基本形式及其调用 4)字符、单元数组、结构体的定义 5)Matlab的流程控制语句第二讲 非线性方程（组）求解与迭代法隋志军 化工学院软件应用教科组 2006-10本章知识要点? 数值计算C 单个非线性方程求解 C 非线性方程组求解 C 迭代法? MATLABC 求解非线性方程（组）的相关函数本章的所要解决的典型问题在945.36kPa（9.33atm）、300.2K时，容器中充以 2mol氮气，试求容器体积。 已知此状态下氮气的P-V-T关系符合范德华方程，其 范德华常数为a＝4.17atm?L/mol2，b＝0.0371L/mol 数学模型：an 2 f (V ) = ( p + 2 )(V ? nb) ? nRT = 0 V非线性方程（组）在化学工程中的作用?多组分混合溶液的沸点、饱和蒸气压计算 ? 流体在管道中阻力计算 ?多组分多平衡级分离操作模拟计算 ?平衡常数法求解化学平衡问题 ?定态操作的全混流反应器的操作分析求解非线性方程的方法? 二分法 ? 不动点迭代 ? 威格斯坦法迭代 ? 牛顿法 ? 割线法迭代法原理不动点迭代： f ( x ) = 0 → x = ? ( x )x 已知： 3 ? x ? 1 = 0的解为1.324，以1.5 为初值，采用以下 两种迭代格式计 算，结果如何？迭代法意义示意图3 x k +1 = x k ? 1 x k +1 = ( x + 1)1 / 3迭代法的收敛性迭代法意义示意图Wegstein法意义示意图牛顿法意义示意图割线法意义示意图多项式函数多项式：P ( x ) = a o x + a1 xnn ?1+ L + a n ?1 x + a nMatlab表达多项式方法： P = [a0 , a1 ,L, a n ?1 , a n ] Matlab多项式函数C C C C C Polyval Polyder Polyfit Conv，Deconv Roots 多项式的值 多项式微分 多项式拟合 多项式乘法和除法 多项式求根简单多项式函数的使用a=[1 11 55 125]; b=[1 1;1 1]; c=polyval(a,b) q=[2,3,5]; aq=conv(a,q) d=poly2sym(aq) e=deconv(aq,a) ％多项式乘法 ％多项式向量表示为符号多项式 ％多项式除法 ％求多项式a在b处的值多项式求根函数roots求解方程：2x 4 ? 5x 3 + 6x 2 ? x + 9 = 0p=[2 -5 6 -1 9]; sol=roots(p) 结果： sol = 1.6024 + 1.4 - 1.2709i -0.3524 + 0.9755i -0.3524 - 0.9755ian 2 f (V ) = ( p + 2 )(V ? nb) ? nRT = 0 VP = 9.33; T = 300.2; n = 2; a = 4.17; b = 0.0371; R = 0.08206 coef=[P, -(P*n*b+n*R*T), a*n^2 , a*b*n^3]; V=roots(coef) 结果： V =5.9 0.1092非线性方程求解函数fzero调用格式： [x,fval,exitflag,output] = fzero(fun,x0,options, p1, p2, ...) 此函数的作用求函数fun在x0附件的零值点，x0是标量。 fval 函数在解x处的值 exitflag 程序结束情况，&0，程序收敛于解； &0，程序没有收敛；＝0，计算达到 了最大次数 output 一个结构体，提供程序运行的信息； output.iterations，迭代次数； output.functions，函数fun的计算次数； output.algorithm，使用的算法 options 选项，可用optimset函数设定选项的新值 fun可以是函数句柄或匿名函数。fzero函数的使用1) sinx在3附近的零点 2) cosx在[1,2]范围内的零点 3) x ? 2 sin x = 03fzero(@sin,3) fzero(@cos,[1,2]) fzero(@(x) x^3-2*sin(x),1) 或者： function Cha2demo1 x=fzero(@fun,1) function y=fun(x) y=x^3-2*sin(x);4) x 3 ? 2 x ? 5 = 0fzero(@(x) x^3-2*x-5,1); roots([1 0 -2 -5])fzero函数初值的选取an 2 f (V ) = ( p + 2 )(V ? nb) ? nRT = 0 VP = 9.33; T = 300.2; n = 2; a = 4.17; b = 0.0371; R = 0.08206; V=fzero(@(V) P*V^3(P*n*b+n*R*T)*V^2+ a*n^2*V a*b*n^3,0)V=0.1092fzero函数初值的选取f (t ) = (sin 2 t )e ? at ? b t的零点以t为自变量，取值范围为-10&t&10，a,b为参数，本例取值分 别为0.1，0.5function Cha2demo2 a=0.1;b=0.5;t=-10:0.01:10; Y=sin(t).^2.*exp(-a*t)-b*abs(t); clf,plot(t,Y,'r');hold on,plot(t,zeros(size(t)),'k'); xlabel('t');ylabel('y(t)'),hold off zoom on n=input('How many zero points are there?'); [tt,yy]=ginput(n);zoom off for i=1:n [t0(i),y(i),exitflag]=fzero(@(t) sin(t)^2*exp(-a*t)-b*abs(t),tt(i)); end disp('The zero points are:') fprintf('%.4f\t',t0) fprintf('\n')fzero函数初值的选取在945.36kPa（9.33atm）、300.2K时，容器中充以2mol氮 气，试求容器体积。 已知此状态下氮气的P-V-T关系符合范德华方程，其范德华常 数为a＝4.17atm?L/mol2，b＝0.0371L/molfunction PVT P = 9.33; % atm T = 300.2; % K n = 2; % mol a = 4.17; b = 0.0371; R = 0.08206; V0 = n*R*T/P [V,fval] = fzero(@PVTeq,V0,[],P,T,n,a,b,R) % -----------------------------------------------------------------function f = PVTeq(V,P,T,n,a,b,R) f = (P + a*n^2/V^2) * (V-n*b) - n*R*T;fzero函数初值对结果的影响function Conv=Cha2demo4 T0 = 450; x0 = [0:0.1:1.0]; n=length(x0); for i=1:n x(i)= fzero(@NonlinEq,x0(i),[],T0); end disp('The conversion could be') fprintf('%.4f\t',x) fprintf('\n') % -----------------------------------------------------------------function f = NonlinEq(x,T0) T = T0 + 250*x; f = 0.25*(1-x)^2*exp(20-10000/T) -结果： The conversion could be 0.8 0.4 0.4 0.3 0.3 1.0951fsolve函数? 与fzero函数只能求解单个方程的根不同，fsolve函数可求 解非线性方程组的解。其算法采用的是最小二乘法。 ? 调用格式： [x,fval,exitflag,output,jacobian] = fsolve(fun,x0,options, p1, p2, ...) ? 输入输入变量的意义同fzero函数。输出变量中的jacobian 为函数fun在x处的Jacobian矩阵。fsolve函数的使用?sin x + y 2 + ln z = 0 ? y 3 ?3x + 2 ? z + 1 = 0 ?x + y + z = 5 ?function Cha2demo6 x0=[1 1 1]; x=fsolve(@fun,x0) function y=fun(x) y(1)=sin(x(1))+x(2)^2+log(x(3))-7; y(2)=3*x(1)+2^x(2)-x(3)^3+1; y(3)=x(1)+x(2)+x(3)-5; 解得结果如下：x＝0.92.0050fsolve函数的应用在铜管内在1 atm下将异丙醇加热到400℃。已知铜是生产丙酮和丙醛的 催化剂，或许还有某些异丙醇异构化为正丙醇。这三种产物的生成可用 如下三个独立反应表示： iC3H7OH(IP) = n C3H7OH(NP) K1 = 0.064 iC3H7OH(IP) = (CH3)CO(AC)+H2 K2 = 0.076 iC3H7OH(IP) = C２H５ＣHO(PR) +H2 K3 = 0.00012 后续加工步骤需要正丙醇，虽然可含丙酮，但丙醛含量不能超过 0.05(mol)%。在上述反应条件下，是否存在违反这种规定的可能性？ 数学模型：各反应的化学平衡方程如下x1 = 0.064 1 ? x1 ? x 2 ? x 3x2 ( x2 + x3 ) = 0.076 (1 ? x1 ? x 2 ? x 3 )(1 + x 2 + x 3 )x3 ( x2 + x3 ) = 0.00012 (1 ? x1 ? x 2 ? x 3 )(1 + x 2 + x 3 )function Cha2demo7 x0 = [0.05 0.2 0.01]; x = fsolve(@EquiC3,x0); CAC=x(3)/sum(x) if CAC&0.05 disp('The AC concentration could not be over 0.05%') else disp('The AC concentration could be over 0.05%') end function f = EquiC3(x) f1 = x(1)-0.064*(1-x(1)-x(2)-x(3)); f2 = x(2)*(x(2)+x(3))-0.076*(1-x(1)-x(2)-x(3))*(1+x(2)+x(3)); f3 = x(3)*(x(2)+x(3))-0.00012*(1-x(1)-x(2)-x(3))*(1+x(2)+x(3)); f = [f1 f2 f3]; 结果：The AC concentration could not be over 0.05%
function Cha2demo5 rho = 961;G = 6.67;L = 10;eps = 5e-6;deltaP = 15e3;k = 1.8;n = 0.64;K = 1.48;n1 = K1 = K*((3*n+1)/(4*n))^n; D1 = (32*K1*L*8^(n1-1)*(4*G/pi/rho)^n1)^(1/(3*n1+1)); D2 = 2*D1; delta = 1e-4; while abs((D2 - D1)/D2) & delta u = (G/rho)/(pi*D1^2/4); Regen = D1^n1*u^(2-n1)*rho/(8^(n1-1)*K1); if Regen &= 2100 f = 16/R end if Regen & 2100 f = FrictFactor(Regen,n1); end Le = k*D1/(4*f); D1 = D2; D2 = (2*f*(L+Le)/(rho*deltaP)*(4*G/pi)^2)^0.2; end fprintf('\t管道直径为D = %.4f %s\n',D2,'m') fprintf('\t摩擦因子为f = %.4f %s',f,'m') % -----------------------------------------------------------------function f = FrictFactor(Regen,n1) f0 = 16/R %4.53 f = fzero(@fFunc,f0,[],Regen,n1) % -----------------------------------------------------------------function F = fFunc(f) F = 1/sqrt(f) - 4.0*log(Regen*f^(1-n1/2))/n1^0.75 + 0.4/n1^1.2;%4.54本讲小结sol=roots(C) [sol,feval,exitflag,output]=fzero(@fun, x0,options,p1,p2,...) [sol,feval,exitflag,output,jacobian]=fso lve(@fun,x0,options,p1,p2,...)第三讲 矩阵操作与线性方程组求解隋志军 化工学院软件应用教科组 2006-10本章知识要点? 数值计算C 线性方程组的直接解法 C 线性方程组的迭代解法? MATLABC 矩阵的生成与操作 C 矩阵运算函数本章的所要解决的典型问题数学模型：∑ ( F (in) ? x (in)) = ∑ ( F (out ) ? x (out )i i线性方程组在化工模型中的作用?多组分体系的物料衡算 ?各种化合物的物理化学性质 ?稳态动力学计算 ?微分方程的差分方程求解病态线性方程组?2 x + 6 y = 8 ? ?2 x + 6.00001 y = 8.00001x=1,y=1?2 x + 6 y = 8 ? ?2 x + 5.99999 y = 8.00002x=10,y=-2?0.2161x + 0.1441 y = 0.1440 ? ?1.2969 x + 0.8648 y = 0.8642x=2,y=-2?0.2161x + 0.1441 y = 0. ? ?1.2969 x + 0.8648 y = 0.x=0.991,y=-0.487线性方程组的求解方法直接求解方法： ? 高斯消元法 ? 高斯－约当消去法 ? 追赶法 迭代解法：C C C C 雅可比（Jacobian）迭代 高斯-赛德尔（Gauss- Siedel）迭代 松弛（SOR）迭代 共轭梯度（CG）迭代在Matlab线性方程组的解法归结为矩阵的除法命令矩阵的生成?直接输入小矩阵 ?利用Array Editor ?创建M文件输入大矩阵 ?利用矩阵操作命令 C cat - 数组连接 C reshape C repmat - 数组变维 - 数组复制与平铺常用矩阵生成函数 工具矩阵zeros ones eye rand randn - 全零阵 - 全一阵 - 单位阵 - 均匀分布随机阵 - 正态分布随机阵B=zeros(n) B=zeros(m,n) B=zeros(size(A))rand(‘state’,n) rand(N) ‘state’ ’twister’ ’seed’常用矩阵生成函数 特殊矩阵compan gallery hadamard hankel hilb invhilb magic pascal rosser toeplitz vander wilkinson - Companion matrix. - Higham test matrices. - Hadamard matrix. - Hankel matrix. - Hilbert matrix. - Inverse Hilbert matrix. - Magic square. - Pascal matrix. - Classic symmetric eigenvalue test problem. - Toeplitz matrix. - Vandermonde matrix. - Wilkinson's eigenvalue test matrix.一维数组的寻访和赋值生成一（1×5）0状态下的均匀分布随机数组，并寻访1)数 组的第三个元素；2) 数组的1，2，5个元素组成的子数组； 3)前三个元素组成的子数组；4)寻访除前三个元素外的其它 元素；5)前三个元素倒排形成的子数组；6)由大于0.5的元 素构成的子数组；7)将元素的第1，4个元素赋值为0。 rand('state',0) x=rand(1,5) a1=x(3) a2=x([1 2 5]) a3=x(1:3) a4=x(3:end) a5=x(3:-1:1) a6=x(find(x&0.5)) x([1 4])=[0 0];x矩阵的生成与寻访A=[1 2 3;4 5 6;7 8 9]; A(5,5)=111 A(:,6)=222 AA=A(:,[1:6,1:6]) AA1=repmat(A,1,2) AA2=repmat(A,2,2) B=ones(2,6);AB_R=[A;B] AB_C=[A,B(:,1:5)'] s=[1 3 6 8 9 11 14 16]; A(s)=0;AA =[ 0 2 0 0 0 ]0 5 0 0 00 8 9 0 00 0 222 0 0 222 0 0 222 0 0 222 0 111 222逻辑1标识逻辑逻辑、关系运算any all find －若向量或数组的列向量任意元素不为0则返回真 －若向量或数组的列向量所有元素不为0则返回真 －寻找非零元素坐标a=magic(5); a (:,3)=zeros(5,1) a 1=all(a(:,1)&10) a 2=all(a&3) a 3=any(a(:,1)&10) a 4=any(a&10)a=reshape([1:9],3,3); a=1./a f1=find(a) f2=find(abs(a)&0.20&abs(a)&0.40)二维数组的操作? 2 6 5 22 2 ? A = ? 1 6 3 8 94? ? ? ?13 4 5 7 21? ? ?求1)数组的第7个元素；2)绝对值最大的元素及其位置；3)所有绝对值大于10的元素 A=[2 6 5 22 2; 1 6 3 8 94; 13 4 5 7 21]; a1=A(7); fprintf('\nThe 7th element is %2d\n',a1) a2=max(max(A)); fprintf('\nThe maximum element of the matrix is %2d\n',a2) [R V]=find(A==a2); fprintf('\nThe maximun element locates at %2dth Row and %2dth colum\n',R,V) L=abs(A)&10; X=A(L)矩阵的特殊操作1)矩阵变维? reshape(X,M,N)2)矩阵变向? ? rot90（A），rot90（A,K） fliplr(A)，flipud(A)，flipdim(X,dim)3)矩阵的抽取? ? diag(X,K) tril(X)，tril(X,K)，triu(X)，triu(X,K)矩阵的基本性质size length ndims numel isempty isequal isequalwithequalnans- Size of array. - Length of vector. - Number of dimensions. - Number of elements. - True for empty array. - True if arrays are numerically equal. - True if arrays are numerically equal.矩阵的基本生成、性质、操作命令归类于 elmat主题，可通过&&help elmat 查看Matlab矩阵函数矩阵分析norm rank det \和/ inv cond lu eig svd －矩阵或向量范数； －矩阵秩； －矩阵的行列式 －线性方程求解 －矩阵求逆 －矩阵条件数 －LU分解 －矩阵的特征值和特征向量 －SVD分解线性方程组求解特征值和特征向量矩阵的分析等较为复杂的函数归类于matfun主 题，可通过&&help matfun查看更为详细的信息左除和右除的比较用以下程序生成条件数很大的方程 randn('state',0);A=gallery('randsvd',100,2e13,2);%产生条件数为2e13的100阶矩阵 x=ones(100,1); %生成真解 b=A*x; %用A和x生成b 比较左除法和求逆法的求解所用的时间和相对残差tic %求逆法 xi=inv(A)*b; ti=toc,eri=norm(x-xi),rei=norm(A*xi-b)/norm(b) tic %左除法 xd=A\b; td=toc,erd=norm(x-xd),rei=norm(A*xd-b)/norm(b) 结果： ti = 0.1950, eri = 0.0539, rei =0.0037 td =0.0156, erd =0.0779, rei =2.本章开始问题的求解数学模型：?0.07 D1 + 0.18 B1 + 0.15 D 2 + 0.24 B 2 = 0.15 × 70 ?0.04 D1 + 0.24 B1 + 0.10 D 2 + 0.65 B 2 = 0.25 × 70 ? ? ?0.054 D1 + 0.42 B1 + 0.54 D 2 + 0.10 B 2 = 0.40 × 70 ?0.035 D1 + 0.16 B1 + 0.21D 2 + 0.01B 2 = 0.20 × 70 ?求解程序：A=[0.07 0.18 0.15 0.24; 0 .04 0.24 0.10 0.65; 0 .54 0.42 0.54 0.10; 0 .35 0.16 0.21 0.01]; B =[0.15*70; 0.25*70; 0.40*70; 0.20*70]; x =A\B反应矩阵的秩与独立反应数独立反应：在一个存在多个反应的复杂反应体系中，某些反应的化学计量方程可以由其它反应的化学计量方程线性组合得到。如果m个同时发生的反应中，若每 一个反应的计量方程都不能有其它反应的计量方程的线性组合得到，则称这m个反 应是相互独立的。反应矩阵：4 NH 3 + 5O2 → 4 NO + 6 H 2 O 4 NH 3 + 3O2 → 2 N 2 + 6 H 2 O 4 NH 3 + 6 NO → 5 N 2 + 6 H 2 O 2 NO + O2 → 2 NO2 2 NO → N 2 + O2 N 2 + 2O2 → 2 NO2?? 4 ? 5 4 6 0 ?? 4 ? 3 0 6 2 ? ?? 4 0 ? 6 6 5 A=? ? 0 ?1 ? 2 0 0 ?0 1 ?2 0 1 ? ? 0 ? 2 0 0 ?1 ?NH3 O2NO H2O N2 NO20? 0? ? 0? ? 2? 0? ? 2? ?b＝rank(A) 解得：b＝3，因此只 有三个反应是独立的线性方程组的迭代解法雅可比迭代：x ( k +1) = D ?1 ( L + U ) x ( k ) + D ?1b高斯-赛德尔迭代：x( k +1)= ( D ? L) Ux?1(k )+ ( D ? L) b?1SOR迭代：x(k+1) =(D?ω )?1[( ?ω)D+ω ]x(k) +ω(D?ω )?1b L 1 U L? 1&ω&2用于加速某收敛的迭代过程，而取0&ω&1用于非收 敛迭代过程使其收敛。 ?最佳松弛因子的选取，一般需在计算过程中搜索寻优线性方程组的迭代解法在偏微分方程差分解法中出现如下典型方程组，试 用分别用雅可比迭代，高斯-赛德尔迭代和SOR法求 解。初值 x ( 0 ) = 10.0, (i = 1,2, L ,10)i? 27 ?? 4 1 ? 1 ?4 1 ? 15 ? ? 1 ?4 1 ? 15 ? O O O M ? ? O O 1 M ? 1 ? 4 ? 15 ? ?? ? ? ? ? ? ? ? ? ?function Cha3demo7 A=[-4 1 0 0 0 0;1 -4 1 0 0 0; 0 1 -4 1 0 0;0 0 1 -4 1 0; 0 0 0 1 -4 1;0 0 0 0 1 -4]; D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1); x0=0*ones(1,6)'; b=[-27 -15 -15 -15 -15 -15]'; %Jacobian Iteration BJ=D\(L+U); gJ=D\b; y=BJ*x0+gJ; n=1; while norm(y-x0)&=1e-6 x0=y; y=BJ*x0+gJ; n=n+1; end disp('The solution of Jacobian iteration are') fprintf('\nn=%3d\n',n) disp('y='),disp(y')%Gauss-Seidel BG=(D-L)\U;gG=(D-L)\b;y=BG*x0+gG; n=1; while norm(y-x0)&=1e-8 x0=y; y=BG*x0+gG; n=n+1; end disp('The solution of Gauss-Seidel iteration are') fprintf('\nn=%3d\n',n) disp('y='),disp(y')%SOR omiga=1.08;BS=(D-omiga*L)\((1omiga)*D+omiga*U);gS=omiga*((D-omiga*L)\b); y=BS*x0+gG; n=1; while norm(y-x0)&=1e-8 x0=y; y=BS*x0+gS; n=n+1; end disp('The solution of SOR iteration are') fprintf('\nn=%3d\n',n) disp('y='),disp(y')第四讲 插值、拟合与数值微分和积分化工学院软件应用教科组 2006-10本章知识要点? 插值方法(interp,spline) ? 拟合方法(polyfit,csaps) ? 数值微分(polyder, fnder) ? 数值积分(quad, quadl, fnint)插值、拟合、数值微分、数值积分 在化工计算中的作用 ?表格式物性数据的内插 ?离散实验数据点的处理 ?状态方程计算流体的焓和熵 ?微分法反应动力学方程拟合 ?等温活塞流反应器的设计计算 ?微观离析反应器的计算插值简介插值的数学问题可以描述为：已知n+1个数对{xi, f(xi)}，其 中i＝0,1…n，（xi互不相同，称之为节点），求取函数 g(xi)=f(xi)。 当{xi, f(xi)}有相当的精确度，但它们的函数关系难以确定或 难以计算时，则可利用这些数据点来构造一个较简单的函 数来近似表达原函数关系。 根据逼近函数的不同，常见的插值方法：?Lagrange多项式插值（线性插值） ?分段插值 ?三次样条插值 ?三角插值 ?有理式插值插值方法的选择已知熔盐在423～473K的密度和粘度如下表所示， 估计450K时的密度和粘度。density~temperature温度 (K) 423 433 443 453 463 473密度 (kg/m3) 59 34粘度 (Pa?S) 177.58 146.51 122.79 104.6 90.26 78.7960 45 30 420430440450460470480viscosity~temperature1801601401201008060 420430440450460470480Matlab的插值（Interpolation）函数插值方法 一维插值 快速一维插值 二维插值 三维插值 N维插值Matlab函数 interp1 interpq interp2 interp3 interpnMatlab函数 插值方法 使用FFT方法的一 interpft 维插值 分段三次Hermite插 pchip 值 spline 三次样条插值一维插值interp1调用格式： yi＝interp1(x,y,xi) 已知数据向量 （x,y），计算并返回在插值向量 xi处的函数值 yi=interp1(x,y,xi, ‘method’) yi=interp1(x,y,xi, ‘method’, ‘extrap’) ‘method’用于指定插值算法，其 值可以是： ‘nearest’――最近插值 ‘linear’――线性插值（默认值） ‘spline’――分段三次样条插值 ‘pchip’――分段三次Hermite插值 ‘cubic’――与‘pchip’相同注意： 向量x为单调。若y为矩阵， 则对y的每一列进行插值 向量xi中有元素不在x的范 围内，则对应yi值为NaN ’extrap’用于指定当向量xi 中有元素不在x的范围内 时，采用’method’所指定的 插值算法进行外插计算与之 对应的yi值一维插值方法比较已知x=[0:10]，y=[0 0.3 0.8 0.4 0.4 0.0];(y＝ sinx)，比较一维线性、线性最近、立方和三次样条插值所得xi＝ 0,0.15,0.30,0.45,…,10处的值yi。如果初始数据点为x＝ 0,2,4,…,10，y＝sinx，以上方法插值效果。x=0:1:10; y =[0 0.3 0.8 -0.4 0.657 0.1 -0.5440]; plot(x,y,'co'),hold on fplot(@sin,[0 10]) xi=0:0.15:10; yi=interp1(x,y,xi); plot(xi,yi,'r+'),text(0.9,'线性插值\rightarrow') yi2=interp1(x,y,xi,'nearst'); plot(xi,yi2,'c*'),text(3.537,0.1374,'\leftarrow最近插值') yi3=interp1(x,y,xi,'cubic'); plot(xi,yi3,'md'),text(2.408,0.8333,'\leftarrow三次插值') yi4=interp1(x,y,xi,'spline'); plot(xi,yi4,'kh'),text(4.62,0.8158,'三次样条插值\rightarrow')初始数据对于插值的影响x=0:2:10;y=sin(x); plot(x,y,'go'),hold on ezplot(@sin,[0 10]) xi=0:0.15:10; yi=interp1(x,y,xi); plot(xi,yi,'r+'),text(0.7,'\leftarrow线性插值') yi2=interp1(x,y,xi,'nearst'); plot(xi,yi2,'c*'),text(6.947,-0.258,'\leftarrow最近插值') yi3=interp1(x,y,xi,‘pchip'); plot(xi,yi3,'md'),text(2.408,0.8333,'\leftarrow三次插值') yi4=interp1(x,y,xi,'spline'); plot(xi,yi4,'kh'),text(1.601,1.138,'\leftarrow三次样条插 值')1.5←三 次 样 条 插 值1←三 次 插 值0.5←线 性 插 值0←最 近 插 值-0.5-1012345678910spline与pchipSpline()的调用格式为：yi=spline(x,y,xi) 此函数等同于yi=interp1(x,y,xi, ‘spline’) pp=spline(x,y) 返回三次样条插值的分段多项式 形式的向量 spline函数可以保证插值函数的三阶导数连续 pchip()的调用格式为： yi=pchip(x,y,xi) 此函数等同于yi=interp1(x,y,xi, ‘pchip’) pp=pchip(x,y) 返回三次样条插值的分段多项式形 式的向量pchip与spline的区别pchip与spline的区别：三次样条在相邻的节点上并不保证单 调性；而Hermite分段三次样条则可保证插值的局部单调性 利用点(x=sin(kπ/6),y=cos(kπ/6)，其中k=[0 1 2 3]来逼近单 位圆的前四分之一圆周。比较pchip与spline的差别。1t=linspace(0,pi/2,4) x=cos(t);y=sin(t); xx=linspace(0,1,40); plot(x,y,'s',xx,[pchip(x,y,xx);spline(x,y,xx)]) grid on, axis equal legend('Orignal data','pchip','spline')0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 Orignal data pchip spline00.20.40.60.81其它样条插值函数－csapecsape()的调用格式为： pp=csape(x,y) pp=csape(x,y,conds) conds可为以下字符串： ‘complete’ 指定端点处一阶导数 ‘second’ 指定端点处二阶导数 ‘periodic’ 左右端点处一、二阶导数相等 ‘not-a-knot’ 第二个和倒数第二个节点处三阶导数连续 ‘variational’ 自然边界条件，即端点二阶导数为零 除csape外，Matlab的样条工具箱还提供了其它样条插值函 数，如csapi，spapi等csape的使用设f(x)为区间[0,3]上的函数，剖分节点为xi=[0 1 2 3]；节点上 的函数值为yi=[0 0.5 2 1.5]，左右端点处的一阶导数值分别 为0.2和-1。试求区间[0,3]上满足上述条件的三次样条插值函 数 程序： x=[0 1 2 3]; y=[0.2 0 0.5 2.0 1.5 -1]; pp=csape(x,y,’complete’)几个有用的函数1. [breaks coff k m n]=unmkpp(pp)命令可以看到样 条函数的具体信息 其中coff是一个矩阵，其第i行是第i个三次多项式 的系数，多项式形式如下：si ( x) = a0 ( x ? xi ?1 )3 + a1 ( x ? xi ?1 ) 2 + a2 ( x ? xi ?1 ) + a 32. fnval或ppval可以计算样条函数在指定点的函数值 3. fnder和fnint可分别计算样条函数的导数和积分 4. fnplt可用于绘制样条函数的图形二维插值：interp2调用格式： zi=interp2(x,y,z,xi,yi,’method’)‘method’算法属性值可以是； ‘nearest’――最近插值 ‘linear’――线性插值（默认） ‘spline’――三次样条插值(spline) ‘cubic’――立方插值二维插值函数的使用假设有一组分度系数的“海底深度测量数据”，由以下一段程序 生成：randn('state',2); x=-5:5;y=-5:5;[X,Y]=meshgrid(x,y); Z=-500+1.2*exp(-((X-1).^2+(Y-2).^2))-0.7*exp(-(exp(X+2).^2+(Y+1).^2)); surf(X,Y,Z)，view(-25,25)试由插值方式绘制海底形状图。xi=linspace(-5,5,50); yi=linspace(-5,5,50); [XI,YI]=meshgrid(xi,yi); ZI=interp2(X,Y,Z,XI,YI,'*cubic'); surf(XI,YI,ZI) view(-25,25)-498.5 -499 -499.5 -500 -500.5 -501 5 0 0 -5 -5-498.5 -499 -499.5 -500 -500.5 -501 5 0 0 -5 -555拟合简介? 拟合和插值的区别在于： 1. 拟合时，所得函数不需要过所有数值点 2. 插值函数不宜外推，拟合函数在某些情况则可以 ? 拟合方法中最常用的是最小二乘曲线拟合 ? 最小二乘法的基本思路是使拟合因变量y在给定点xi上 使残差平方和最小 ? 用于拟合的函数可以是多样的，本讲只介绍多项式函 数拟合和样条函数拟合最小二乘多项式拟合：polyfit调用格式如下： p=polyfit(x,y,n) [p,s]=polyfit(x,y,n) 说明： 输入参数：(x,y)为已知数据向量，n为多项式阶数；输出 参数p为拟合生成的多项式的系数向量（长度为n+1），s为 结构参数，供函数polyval()调用以获得误差估计值。 函数polyval()常常与polyfit()联合使用，其调用格式为： y=polyval(p,x)。 [y,deltay]=polyval(p,xi,s), 返回xi处的拟合函数值，deltay 是50%置信区间的y的误差。多项式次数对拟合效果的影响已知下表数据，用polyfit进行多项式拟合x y 0.5 1.75 1.0 2.45 1.5 3.81 2.0 4.80 2.5 8.00 3.0 8.60x=0.5:0.5:3;y=[1.75 2.45 3.81 4.80 8.00 8.60]; xi=0.5:0.05:3; [a2,S2]=polyfit(x,y,2);[y2,delta2]=polyval(a2,xi,S2); plot(x,y,'ro',xi,y2,'m-'),text(1.04,4.67,'二次多项式拟合\leftarrow') hold on, [a4,S4]=polyfit(x,y,4);[y4,delta4]=polyval(a4,xi,S4); plot(x,y,'ro',xi,y4,'b-'),hold on,text(1.563,6.775,'四次多项式拟合\leftarrow') [a7,S7]=polyfit(x,y,7);[y7,delta7]=polyval(a7,xi,S7); plot(x,y,'ro',xi,y7,'k-'),hold on,text(1.931,9.532,'七次多项式拟合\leftarrow'),hold off最小二乘法拟合生成样条曲线函数名 csaps() spap2() spaps() 曲线类型 三次样条曲线 B样条曲线 B样条曲线 拟合准则 最小二乘法 最小二乘法 最小二乘法是否平滑处理是 否 是? 拟合得到曲线函数sp以后，可利用fnval()计算任 意自变量下的函数值。函数csaps()的用法功能：平滑生成三次样条函数，即对于数据(xi,yi)，所求的三次样条 函数y=f(x)满足 ：min p ? ∑ w i ( y i ? f ( x i )) 2 + (1 ? p ) ∫ λ (t ) D 2 f t 2 dti()调用格式：sp=csaps(x,y,p) ys=csaps(x,y,p,xx,w) 输入参数：x,y 要处理的离散数据(xi,yi) p 平滑参数，取值区间为[0,1]。当p=0时，相当于最小二乘直线拟 合；当p=1时，相当于“自然的”三次样条拟合 xx 用于指定在给定点xx上计算其三次样条函数值(ys) w 权值（权重），默认为1 输出参数：sp 拟合得到的样条函数 ys 在给定点xx上的三次样条函数值函数spap2 ()的用法功能：用最小二乘法拟合生成B样条曲线，即对于数据(xi,yi)，所求k次样条函数y=f(x)满足 ：min ∑ wi y i ? f ( xi )i()2调用格式：sp=spap2(knots,k,x,y) sp=spap2(knots,k,x,y,w) knots k x,y w sp 节点序数(knot sequence) 样条函数的阶次，一般取k=3，有时取k=4 要处理的离散数据(xi,yi) 权值（权重），默认为1 拟合得到的样条函数输入参数：输出参数：函数spaps ()的用法功能：平滑生成B样条函数min ∑ wi y i ? f ( xi )i()2调用格式：sp=spaps(x,y,tol) [sp,ys]=spaps(x,y,tol,m,w) x,y tol m w sp ys 要处理的离散数据(xi,yi) 光滑时的允许精度 默认值是2，即平滑生成三次B样条曲线 权值（权重），默认为1 拟合得到的样条函数 在x上经平滑处理的B样条函数值输入参数：输出参数：函数csaps()的用法采用样条拟合函数csaps进行拟合，比较p的取值对于拟合 效果的影响x y 0.5 1.75 1.0 2.45 1.5 3.81 2.0 4.80 2.5 8.00 3.0 8.60x=0.5:0.5:3; y=[1.75 2.45 3.81 4.80 8.00 8.60]; xi=0.5:0.05:3; plot(x,y,'ro') hold on, C=[0 0 0;1 0 0;0 0 1; 0 1 1; 0.5 0.5 0.5]; j=1; for i=0:0.25:1.0 Cj=C(j,:); ys(j,:)=csaps(x,y,i,xi); plot(xi,ys(j,:),'color',Cj),pause j=j+1; end对于微小的数据扰动，多项 式拟合通常比样条函数更为 敏感数值微分? 对于列表型函数往往需要用数值方法计算函数的微分 ? 数值微分的基本方法 1. 差分 2. 利用插值（拟合）多项式求微分 3. 利用三次样条插值（拟合）函数求微分 ? 数值微分可以放大误差，应谨慎使用Matlab数值微分实现方法有限差分法：用差分函数diff()近似计算导数 多项式拟合方法：polyfit () polyder polyval () 离散数据 ?? ?→ 多项式拟合函数 ?? ?() → 导函数pp ?? ? → 在xi的导数值 ? ?三次样条插值方法fnder () fnval () 离散数据 ?spline () → 样条插值函数cs ???→ 导函数pp ???→ 在xi的导数值 ??样条拟合方法csaps()或spap2 a ()或spaps() fnder () fnval () ? 离散数据 ??????? → 样条拟合函数sp ???→导函数pp ???→ 在xi的导数值函数diff对于向量X，diff(X)表示了[X(2)-X(1) X(3)-X(2) ... X(n)-X(n-1)]. 对于矩阵X，diff(X)表示了[X(2:n,:) - X(1:n-1,:)] diff(x,n,dim)得到矩阵x在dim维上的n阶差值&& diff((1:10).^2) x=[1 3 8; 2 4 6] diff(x,1,1) diff(x,1,2) [1 3 5 7 9 11 13 15 17 19] diff(x,2,2) diff(x,3,2)1 1 -2[2 5;2 2][3 0]Empty matrix: 2-by-0利用diff函数求sin(x) 在[0,10]上的导数 值fplot(@cos,[1 10]) hold on h=1;x=1:h:10; hh=0.01;xx=1:hh:10; diffy=diff(sin(x))/h; diffyy=diff(sin(xx))/ plot([1:h:10-h],diffy,'r:',[1:hh:10-hh],diffyy,'k-o')函数gradientFX = gradient(F) [FX,FY] = gradient(F) [Fx,Fy,Fz] = gradient(F) [...] = gradient(F,h) [...] = gradient(F,h1,h2,...) 其中 1. Fx=dF/dx，FY=dF/dy，Fz=dF/dz；当F为向量时，dF=gradient(F)是一维 梯度 2. h1, h2，……是步长，缺省值为1 3. 如果h1，h2，……为向量,其长度必须匹配F的维数例题某一液相反应浓度随时间变化的实验数据如下表：t/minC(g/L)0 5.190.2 3.770.6 2.301.0 1.572.0 0.85.010.00.25 0.094试t=0.1,0.4min时的反应速度x=[0 0.2 0.6 1 2 5 10]; C=[5.19 3.77 2.30 1.57 0.8 0.25 0.094]; plot(x,C,'bo'),hold on pp=csaps(x,C);fnplt(pp) dC=fnder(pp); dC1=ppval(dC,0.1) dC2=fnval(dC,0.4) fprintf('The reaction rate at 0.1min and 0.4min are%5.4f,%5.4f... (g/(L min)respectively\n',dC1,dC2)数值积分? 数值积分在数值计算中有着重要作用,许多数值计算问题 可以转化为数值积分问题,如常微分方程初值问题等 ? 通常可用逼近多项式Pn(x)来代替被积函数f(x)，计算积 分 ∫ P ( x)dxb a n? 构造数值积分的方法很多，主要有Newton-Cotes系列 数值积分法、Gauss积分法和Romberg积分法等数值积分MATLAB函数 quad quadl trapz cumtrapz cumsum fnint 公式 自适应Simpson求积公式（低阶） 自适应Lobatto求积公式；精度高，最常 用 梯形求积公式；速度快，精度差 梯形法求一个区间上的积分曲线 等宽距法求一个区间上的积分曲线，精 度很差 利用样条函数求不定积分；与spline， ppval配合使用，主要应用于表格“函数” 积分梯形法数值积分：trapz()调用格式： z=trapz(y)用梯形求积方法计算y的积分近似值。 对于向量y，trapz(y)返回y的积分； 对于矩阵y，trapz(y)返回一行向量，向量中的各元素为矩 阵y的对应列 向量的积分值；自适应Simpson法数值积分：quad()调用格式： q=quad(@fun,a,b) q=quad(@fun,a,b,tol) q=quad(@fun,a,b,tol,trace,p1,p2,……) fun 被积函数。在定义fun时，被积函数表达 式必须是向量形式，即表达式必须使用点 运算符(.*、./和.^)以支持向量 a,b 即积分限[a,b] tol 绝对误差限，默认值为1.e-6 p1,p2…… 直接传递给函数fun的已知参数 q 积分结果输入参数：输出参数：自适应Lobatto法数值积分：quadl() 调用格式同quad不同积分函数的比较求积分： ∫03πe ?0.5t sin(t + π / 6)dt比较cumsum， trapz，quad， quadl的积分精 度，该积分的精确 解为 0.function Cha4demo5 format long d=pi/1000; t=0:d:3* nt=length(t); y=fun(t); sc=cumsum(y)*d; scf=sc(nt) z=trapz(y)*d qd=quad(@fun,0,3*pi,[],1) qd2=quadl(@fun,0,3*pi,[],1) EV=0.; err(1)=abs(scf-EV); err(2)=abs(z-EV)*10; err(3)=abs(qd-EV)*10; err(4)=abs(qd2-EV)*10; bar(err),title('不同积分方法比较'), colormap(summer) text(0.7,7e-4,'矩形法','fontsize',20) text(1.5,5e-4,'梯形法','fontsize',20) text(2.4,3e-4,'Simpson法','fontsize',20) text(3.4,1e-4,'Lobatto法','fontsize',20) %-----------------------------------------------------------function y=fun(t) y=exp(-0.5*t).*sin(t+pi/6)广义积分1. 奇点积分 2. 无穷积分∫10dx x (exp(x) + 1)∫+∞0exp(sin x ? x 2 / 100)dxfun=inline('1./(sqrt(x).*(exp(x)+1))'); quadl(fun,0,1) quadl(fun,eps,1)可先选取一个有限的积分区间，如 [0,100]计算；在选择一个较大的积 分区间，如[0 200]计算，如两次计 算结果的差满足一定的精度要求， 则可认为此值即为无穷积分的值多重数值积分二重积分函数：SS=dblquad(fun,xmin,xmax,ymin,ymax,tol,method)三重积分函数：SSS=triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method)说明: Matlab只能处理积分限为常数的多重积分，对内积分上下限为 外积分变量的积分问题,需自行编写函数求解。样条函数在数值积分与微分中的应用已知x=(0:0.1:1)*2*y=sin(x) ，求其积分和导数x=(0:0.1:1)*2*y=sin(x); pp=spline(x,y);int_pp=fnint(pp); der_pp=fnder(pp);% 在基础区间上，计算三个样条函数与理论值的最大误差 xx=(0:0.01:1)*2* err_yy=max(abs(ppval(pp,xx)-sin(xx))) err_int=max(abs(ppval(int_pp,xx)-(1-cos(xx)))) err_der=max(abs(ppval(der_pp,xx)-cos(xx))) % 计算y(x)在区间[1,2]上的定积分 DefiniteIntegral.byTheory=(1-cos(2))-(1-cos(1)); % 计算dy(3)/dx Derivative.bySpline=fnval(der_pp,3); Derivative.byTheory=cos(3); Derivative.byDiference=(sin(3.01)-sin(3))/0.01; %前向差分近似 DefiniteIntegral,Derivative fnplt(pp,'b-');hold on fnplt(int_pp,'m:'),fnplt(der_pp,'r--');hold off legend('y(x)','S(x)','dy/dx')ppval(int_pp,[1,2])给出一个行向 量，其中第一个元素是原函数在 0～1的定积分，第2个元素是原 函数在0～2的定积分，所以在 1～2的定积分应是第2个元素减 第1个元素的值DefiniteIntegral.bySpline=ppval(int_pp,[1,2])*[-1;1];反应器停留时间分布的混合特性在t＝0的时刻，在一容器入口处突然向流进容器的流体脉冲注 入一定量的示踪剂，同时在容器出口处测量流出物料中示踪剂 浓度随时间的变化，实验数据如下表：t/sC(kmol/m3)×1030 020 040 060 080 0.4100 5.5120140160 1.7180 0.1200 016.2 11.1试计算流体在容器中的平均停留时间以及扩散准数。 数学模型：平均停留时间 方差： 扩散特征数为：∞tm∫ Ctdt = ∫ Cdt0 ∞ 0σ2∫ Ct dt ? t = ∫ Cdt2 0 ∞ 0∞2 m? DL ? σ 2 ? ?= 2 ? uL ? 2t mfunction Cha4demo7 t = [0:20:200];C = [0, 0, 0, 0, 0.4, 5.5, 16.2, 11.1, 1.7, 0.1, 0];plot(t,C,'o') sp1= spline(t,C); hold on fnplt(sp1,'b-'); xlabel(‘Time (s)’)，ylabel('C (kmol/m^3)??10^3') hold off % By spaps(): figure，plot(t,C,'o') sp2= spaps(t,C,1); hold on fnplt(sp2,'b-'); xlabel(‘Time (s)’)，ylabel('C (kmol/m^3)??10^3') hold off sp=input('Which function should be used to integrate?'); % Integration t0 = 0;tf = t(end); IC = quadl(@Func,t0,tf,[],[],sp); ICt = quadl(@Func1,t0,tf,[],[],sp); ICt2 = quadl(@Func2,t0,tf,[],[],sp); tm1 = ICt/IC ss1 = ICt2/IC - tm1^2 DL2uL1 = ss1/tm1^2/2 % -----------------------------------------------------------------function y = Func(x,sp) % f= C y = fnval(sp,x); % -----------------------------------------------------------------function y = Func1(x,sp) % f= Ct y = fnval(sp,x).*x; % -----------------------------------------------------------------function y = Func2(x,sp) % f= Ct^2 y = fnval(sp,x).*x.^2;微分法进行动力学数据分析A ， 反应物A在一等温间歇反应器中发生反应为：→ 产物 测量得到反应器中不同时间下反应物A的浓度CA如 下表所示。试根据表中数据确定其反应速率方程。t/s CA(mol/L) 0 10 20 8 40 6 60 5 120 3 180 2 300 1n 数学模型： ? rA = kC Aln(?rA ) = n ln C A + ln kfunction Cha4demo8 % 动力学数据 t = [0 20 40 60 120 180 300];CA = [10 8 6 5 3 2 1]; % 用最小二乘样条拟合法计算微分dCA/dt knots = 3;K = 3; % 三次B样条 sp = spap2(knots,K,t,CA); sp = spap2(newknt(sp),K,t,CA); pp = fnder(sp); % 计算B样条函数的导函数 dCAdt = fnval(pp,t); % 计算t处的导函数值 % 绘制图形 ti = linspace(t(1),t(end),200); CAi = fnval(sp,ti); plot(t,CA,‘ro’,ti,CAi,‘b-’)，xlabel(‘t’)，ylabel('C_A') figure fnplt(pp); % dCAdti = fnval(pp,ti) % plot(ti,dCAdti,'-') xlabel('t') ylabel('dC/dt') % 线性拟合 rA = dCA y = log(-rA); x = log(CA); p = polyfit(x,y,1); k = exp(p(2)) n = p(1)第五讲 常微分方程数值解化工学院软件应用教科组 2006-10本章知识要点? 数值计算C 常微分方程初值问题 C 常微分方程边值问题? MATLABC 微分方程求解常微分方程的相关函数? ode45 ode23 ? bvp4c微分方程在化工模型中的应用 ?间歇反应器的计算 ?活塞流反应器的计算 ?全混流反应器的动态模拟 ?定态一维热传导问题 ?逆流壁冷式固定床反应器一维模型 ?固定床反应器的分散模型Matlab常微分方程求解问题分类初值问题： ? 定解附加条件在自变量 的一端 ? 一般形式为：? y ' = f ( x, y ) ? ? y (a) = y 0边值问题： 在自变量两端均给定附加 条件? y ' = f ( x, y ) 一般形式：? y(a) = y1 , y(b) = y2 ?? 初值问题的数值解法一 般采用步进法，如 Runge-Kutta法边值问题可能有解、也可 能无解，可能有唯一解、 也可能有无数解 边值问题有3种基本解法 ? 迭加法 ? 打靶法 ? 松弛法Matlab求解常微分方程初值问题方法 将待求解转化为标准形式，并“翻译” 成Matlab可以理解的语言，即编写 odefile文件 选择合适的解算指令求解问题根据求解问题的要求，设置解算指令 的调用格式Matlab求解初值问题函数指令 ode23 含义 普通2-3阶法解ODE 指令 odefile 含义 ODE文件格式ode45 ode113 解 算 ode23t ode15s ode23s ode23tb普通4-5阶法解ODE 普通变阶法解ODE 解适度刚性ODE 变阶法解刚性ODE 低阶法解刚性ODE 低阶法解刚性ODE 输 出 选 项odeset odeget odeplot odephas2 odephas3 odeprint创建、更改ODE选项的 设置 读取ODE选项的设置 ODE的输出时间序列图 ODE的二维相平面图 ODE的三维相空间图 在Matlab指令窗显示结 果odefile所谓的odefile实际上是一个Matlab函数文件，一般作 为整个求解程序的一个子函数，表示ode求解问题 Matlab提供了odefile的模板，采用type odefile命令显 示其详细内容，然后将其复制到脚本编辑窗口，在合 适的位置填入所需内容 一般而言，对于程序通用性要求不高的场合，只需将 原有模型写成标准形式，然后“翻译”成Matlab语言即可odefile的编写规定1. ode文件的最简单格式必须有一个自变量t和 函数y作为输入变量，一个y的导函数作为输 出变量。其中自变量t不论在ode文件中是否 使用都必须作为第一输入变量，y则必须作为 第二输入变量，位置不能颠倒。 2. 可以向ode文件中传递参数，数目不受限制odefile的编写求解初值问题：2x ? y' = y ? ? y ? ? y (0) = 1 ?(0 ≤ x ≤ 1)? y ' = f ( x, y ) ? ? y (a) = y 0ode输入函数 function f=fun(x,y) f=y-2*x/y; 输出变量为因变 量导数的表达式自变量在前，因变 量在后? y' = y + y 2 初值问题：? ? y ( 0) = 1(0 ≤ x ≤ 1)function f=fun(x,y) f=y＋y^2;常微分方程组odefile的编写? y1 ' = 0.04(1 ? y1 ) ? (1 ? y 2 ) y1 + 0.0001(1 ? y 2 ) 2 ? 4 2 ? y 2 ' = ?10 y1 '+3000(1 ? y 2 ) ? y (0) = 0, y (0) = 1,0 ≤ x ≤ 100 2 ? 1常微分方程组与单个常微分方程求解方法相同， 只需在编写odefile时将整个方程组作为一个向量 输出。function f=fun(x,y) dy1dx = 0.04*(1-y(1))-(1-y(2)).*y(1)+0.0001*(1-y(2)).^2; dy2dx = -1e4*dy1dx + 3000*(1-y(2)).^2; f = [dy1 dy2dx];高阶微分方程odefile的编写求解： y&+ a (t )( y ' ) 2 + b(t ) y = e t cos 2πta (t ) = ?e ? t + cos 2πte ?2t , b(t ) = cos(2πt )y(0)=0，y'(0)=1，本例的难度： 方程系数非线性 方程高阶，非标准形式 方程变形：令y1＝y；y2＝y’ 则原方程等价于：function f=fun(t,y) a=-exp(-t)+cos(2*pi*t)*exp(-2*t); b=cos(2*pi*t); f=[y(2) -a*y(2)^2-b*y(1)+exp(t)*b];? y1' = y 2 ? ? ' 2 ? y 2 = ? ay 2 ? by1 + e t cos 2πt ?可在odefile中定义解算指令的使用方法调用格式：1. 2. 3. 4. [T,Y]=ode45(@fun, TSPAN,Y0) [T,Y]=ode45(@fun, TSPAN,Y0,options) [T,Y]= ode45(@fun, TSPAN,Y0,options,P1,P2,…) [T,Y,TE,YE,IE]= ode45(@fun, TSPAN,Y0,options,P1,P2,…)说明：输出变量T为返回时间列向量；解矩阵Y的每一行对应于T的一个元素，列数与求解 变量数相等。 @fun为函数句柄，为根据待求解的ODE方程所编写的ode文件（odefile）； TSPAN＝[T0 TFINAL]是微分系统y'＝F(t,y)的积分区间；Y0为初始条件 options用于设置一些可选的参数值，缺省时，相对于第一种调用格式。options中 可以设置的参数参见odeset P1，P2，…的作用是传递附加参数P1，P2，…到ode文件。当options缺省时，应 在相应位置保留[]，以便正确传递参数。常微分方程初值问题解算指令比较解算指令 ode45 ode23 ode113 ode15s ode23s ode23t ode23tb 算法 四五阶Runge-Kutta法 二三阶Runge-Kutta法 可变阶Adams-Bashforth-Moulton法 基于数值差分的可变阶方法（BDFs， Gear） 二阶改进的Rosenbrock法 使用梯形规则 TR-BDF2（隐式Runge-Kutta法） 低～中 低 适中 低 精度 较高 低ode解算指令的选择(1)1.根据常微分方程要求的求解精度与速度要求2x ? y' = y ? ? y ? ? y ( 0) = 1 ?求解初值问题：(0 ≤ x ≤ 1)比较ode45和ode23的求解精度和速度ode45和ode23的比较function Cha5demo1 %Comparison of results obtained from ode45 and ode23 solver format long y0=1; tic,[x1,y1]=ode45(@fun,[0,1],y0);t_ode45=toc tic,[x2,y2]=ode23(@fun,[0,1],y0);t_ode23=toc plot(x1,y1,'b-',x2,y2,'m-.'), xlabel('x'),ylabel('y'), legend('ODE45','ODE23','location','Northwest') disp('Comparative Results at x=1:'); fprintf('\nODE45\t\t\t y=%.8f\nODE23\t\t\t y=%.8f\nPrecisive Result ... =%.8f\n',y1(end),y2(end),1.7320508) %-----------------------------------------------------------------------function f=fun(x,y) f=y-2*x/y;ode解算指令的选择(2)2.根据常微分方程组是否为刚性方程? y ' = Ay + b( x) ? ? y ( x0 ) = y0如果系数矩阵A的特征值连乘积小于零， 且绝对值最大和最小的特征值之比（刚性 比）很大，则称此类方程为刚性方程 刚性比：100/0.01＝10000? y1' = ?0.01y1 ? 99.99 y2 ? ' ? y2 = ?100 y2 ? y (0) = 2, y (0) = 1 2 ? 1? 刚性方程在化学工程和自动控制领域的模型中比较常 见。ode解算指令的选择(2)? 刚性方程的物理意义：常微分方程组所描述的物理化学变 化过程中包含了多个子过程，其变化速度相差非常大的数 量级，于是常微分方程组含有快变和慢变分量。 ? 常微分方程组数值积分的稳定步长受模值最大的特征值控 制，即受快变量分量约束，特征值大则允许步长小；而过 程趋于稳定的时间又由模值最小的特征值控制，特征值小 则积分到稳定的时间则长。 ? Matlab提供了不同种类的刚性方程求解指令：ode15s ode23s ode23t ode23tb，可根据实际情况选用刚性常微分方程组求解function Cha5demo3 figure ode23s(@fun,[0,100],[0;1]) hold on, ode45(@fun,[0,100],[0;1]) %-------------------------------------------------------------------------function f=fun(x,y) dy1dx = 0.04*(1-y(1))-(1-y(2)).*y(1)+0.0001*(1-y(2)).^2; dy2dx = -1e4*dy1dx + 3000*(1-y(2)).^2; f = [dy1 dy2dx];解算指令的输出控制1. [T,Y]= ode45(@fun, TSPAN,Y0,options,P1,P2,…) 2. sol=ode45(@fun,TSPAN,Y0)将解输出指结构体sol中；SXINT = deval(SOL,XINT)用于计算解在XINT处的值，XINT必须位于区间[SOL.x(1) SOL.x(end)]内。 3. 在无输出变量时，将调用默认的odeplot输出解的图形。 4. 除了以odeplot形式输出外，还可以以odephas2，和odephas3的形式输出解 向量的二维和三维相平面图。 5. 采用以下语句options=odeset(‘outputfcn’,‘odephas2’)可以将输出方法改变为 平面输出 6. odeprint输出求解过程每一步的解解算指令的options选项1. RelTol－相对误差，它应用于解向量的所有分量。在每一步积分过程中，第i 个分量误差e(i)满足：e(i)&=max(RelTol*abs(y(i),AbsTol(i))。 2. AbsTol－绝对误差，若是实数，则应用于解向量的所有分量，若是向量，则 它的每一个元素应用于对应位置解向量元素。 3. OutputFcn－可调用的输出函数名。每一步计算完后，这个函数将被调用输 出结果，可以选择的值为：odeplot，odephas2，odephas3，odeprint。 4. OutputSel－输出序列选择。指定解向量的哪个分量被传递给OutputFcn。 5. MaxSetp－步长上界，缺省值为求解区间的1/10。 6. InitialStep－初始步长，缺省时自动设置。 7. 采用odeset改变原有选项的值间歇反应器中串联－平行复杂反应系统在间歇反应器中进行液相反应制备产物B，反应网络如图5-1所示。反 应可在180~260℃的温度范围内进行，反应物X大量过剩，而C, D和 E为副产物。各反应均为一级动力学关系：r＝－kC，式中 k = k e? Ea RT 0已知参数：k01＝5.7，k02＝3.9，k03＝ 1.6，k04＝6.264×108，Ea1=124670，Ea2=150386， Ea3=77954，Ea4=111528。初始浓度：CA=1kmol/m3，其余物质浓度为0。 已知是产物B收率最大的最优反应温度为224.6℃ 试计算1)在最优反应温度下各组分浓度随时间的动态变化；2)最优反应时 间；3)输出产物D对反应物浓度A的关系图。间歇反应器中串联－平行复杂反应系统数学模型dC A = ?(k1 + k 2 )C A dtdC B = k1C A ? k 3C B dtdCC = k 2 C A ? k 4 CC dt dC D = k 3C B ? k 5 C D dt dC E = k 4 CC + k 5C D dt间歇反应器中串联－平行复杂反应系统function Cha5demo4 T = 224.6 + 273.15; R = 8.31434; k0 = [5.7.9.6.264E+8]; Ea = [386 ]; % Initial concentration:C0(i), kmol/m^3 C0 = [1 0 0 0 0]; tspan = [0 1e4];opt=odeset(‘reltol’,1e-4,’outputfcn’,’odephas2’,’outputsel’,[1;4])[t,C] = ode45(@MassEquations, tspan, C0,opt,k0,Ea,R,T) % plot plot(t,C(:,1),'r-',t,C(:,2),'k:',t,C(:,3),'b-.',t,C(:,4),'k--'); xlabel('Time (s)'); ylabel('Concentration (kmol/m^3)'); legend('A','B','C','D') CBmax = max(C(:,2)); % CBmax yBmax = CBmax/C0(1) % yBmax: index = find(C(:,2)==CBmax); t_opt = t(index) % t_opt: the optimum batch time, s % -----------------------------------------------------------------function dCdt = MassEquations(t,C,k0,Ea,R,T) % Reaction rate constants, 1/s k = k0.*exp(-Ea/(R*T)); k(5) = 2.16667E-04; % Reaction rates, kmoles/m3 s rA = -(k(1)+k(2))*C(1); rB = k(1)*C(1)-k(3)*C(2);rC = k(2)*C(1)-k(4)*C(3); rD = k(3)*C(2)-k(5)*C(4); rE = k(4)*C(3)+k(5)*C(4); % Mass balances dCdt = [rA; rB; rC; rD; rE];Matlab求解边值问题方法：bvp4c函数1. 把待解的问题转化为标准边值问题? y ' = f ( x, y ) ? ? g ( y (a), y (b)) = 02. 因为边值问题可以多解，所以需要为期望解指定一个初 始猜测解。该猜测解网（Mesh）包括区间[a, b]内的一组 网点（Mesh points）和网点上的解S(x) 3. 根据原微分方程构造残差函数r ( x) = S ' ( x) ? f ( x, S ( x))4. 利用原微分方程和边界条件，借助迭代不断产生新的 S(x)，使残差不断减小，从而获得满足精度要求的解Matlab求解边值问题的基本指令solinit＝bvpinit（x,v,parameters） 生成bvp4c调用指令所必须的“解猜测网” sol＝bvp4c（odefun，bcfun，solinit，options，p1，p2，…） 给出微分方程边值问题的近似解 sxint＝deval（sol，xint） 计算微分方程积分区间内任何一点的解值初始解生成函数：bvpinit()solinit＝bvpinit（x,v,parameters）x指定边界区间[a,b]上的初始网络，通常是等距排列的（1×M）一维数组。 注意：使x(1)=a，x(end)=b；格点要单调排列。 v是对解的初始猜测 solinit（可以取别的任意名）是“解猜测网（Mesh）”。它是一个结构体， 带如下两个域： solinit.x是表示初始网格有序节点的（1×M）一维数组，并且 solinit.x(1)一定是a，solinit.x(end)一定是b。M不宜取得太大，10数 量级左右即可。 solinit.y是表示网点上微分方程解的猜测值的（N×M）二维数组。 solinit.y(:,i)表示节点solinit.x(i)处的解的猜测值。边值问题求解指令：bvp4c ()sol＝bvp4c（odefun，bcfun，solinit，options，p1，p2，…） 输入参数：odefun是计算导数的M函数文件。该函数的基本形式为：dydx＝ odefun(x,y,parameters,p1,p2，…)，在此，自变量x是标量，y，dydx 是列向量。其基本形式与初值问题的odefile形式相同 bcfun是计算边界条件下残数的M函数文件。其基本形式为：res＝ bcfun(ya,yb,parameters,p1,p2,……)，文件输入宗量ya,yb是边界条件 列向量,分别代表y在a和b处的值。res是边界条件满足时的残数列向量。 注意：例如odefun函数的输入宗量中包含若干“未知”和“已知”参数，那 么不管在边界条件计算中是否用到，它们都应作为bcfun的输入宗量。 输入宗量options是用来改变bvp4c算法的控制参数的。在最基本用法 中，它可以缺省，此时一般可以获得比较满意的边值问题解。如需更改 可采用bvpset函数，使用方法同odeset函数。 输入宗量p1，p2等表示希望向被解微分方程传递的已知参数。如果无 须向微分方程传递参数，它们可以缺省。边值问题求解指令：bvp4c ()输出参数： 输出变量sol是一个结构体 sol.x是指令bvp4c所采用的网格节点； sol.y是y(x)在sol.x网点上的近似解值； sol.yp是y'(x)在sol.x网点上的近似解值； sol.parameters是微分方程所包含的未知参数的近似解值。当 被解微分方程包含未知参数时，该域存在。边值问题的求解? y&? y = 10 求解两点边值问题：? ? y (0) = y (1) = 0令： y1 = y, y 2 = y1 '原方程组等价于以下标准形式 的方程组：? y1 ' = y 2 ? ? y 2 ' = y1 + 10边界条件为：y1 (0) = 0, y1 (1) = 0function Cha5demo5 solinit=bvpinit(linspace(0,1,10),[1 0]); sol=bvp4c(@ODEfun,@BCfun,solinit); x=[0:0.05:0.5]; y=deval(sol,x); yP=[0 -0..... -0...]; xP=[0:0.1:0.5]; plot(xP,yP,'o',x,y(1,:),'r-'), legend('Analytical Solution','Numerical Solution') % -----------------------------------------------------------------function dydx=ODEfun(x,y) dydx=[y(2);y(1)+10]; % -----------------------------------------------------------------function bc=BCfun(ya,yb) bc=[ya(1);yb(1)];边值问题的求解求解：? y&?(1 + x 2 ) y = 1 ? ? y (0) = 1, y (1) = 3function Cha5demo6 solinit = bvpinit(linspace(0,1,10),[0 1]); sol = bvp4c(@ODEfun,@BCfun,solinit); x = [0:0.1:1]; y = deval(sol,x); yP=[1 1.5 1.4... 1.7 2.4 2.6214 3]; xP=[0:0.