Matlab与化学化工计算

分段三次 Hermite插值
5
6
7
-1
0
1
2
3
4
5
6
7
例7-5

用函数y=ex生成以下离散数据，使用Matlab的 不同插值方法计算x=[2.55 2.63 2.77 2.86]处的 函数值，并与真实值进行比较
x 2.5 2.6 2.7 2.8 2.9
y
12.1825
13.4637
14.8797
16.4446
2 线性方程组求解
2.1 线性方程组的一般形式

在应用中，常常把线性方程组
a11 x1 a12 x2 a1n xn b1 a x a x a x b 21 1 22 2 2n n 2 as1 x1 as 2 x2 asn xn bs
计算机在化学化工中的应用
七 Matlab与化学化工计算
本节要点



本章背景 Matlab基础 方程组求解 数据插值 作业
问题的提出


MATLAB语言与其它语言的关系仿佛和 C语言与汇编语言的关系一样 计算机语言的发展
管理、可视化 智能化 解析运算
数值运算

标志着计算机语言向“智能化”方向发 展，被称为第四代编程语言
Matlab 的发展
1984年，MATLAB第1版(DOS版)
1992年，MATLAB 4.0版 1994年，MATLAB 4.2版 1997年，MATLAB 5.0版 1999年，MATLAB 5.3版
告别DOS版 1993年 MathWorks 公司从加拿大 滑铁卢大学购 的 的使 得5.0 Maple MATLAB 拥有 用权，推出了 更丰富的数据类 符号计算工具 型和结构、更友 包 善的面向对象、 更加快速精良的 图形可视、更广 博的数学和数据 分析资源、更多 的应用开发工具
y2 y1 y=?
y0

x0
x1
x2
x
y2 y0 y1 y0 y1 y0 x2 x0 x1 x0 y y0 ( x x0 ) ( x x0 )(x x1 ) x1 x0 x2 x1
3.3.1 使用Matlab进行数据插值

一维插值只有一个自变量的插值 Matlab提供的一维插值函数是interp1
“%”是 Matlab的注释 符，%后的语 句作为注释处 理
解：>> a=[7 0 28; 0 28 1; 28 0 196]; >> b=[1 -39 -7]'; %b为列向量，故输入行向量后转置 >> rank(a) %计算系数矩阵A的秩 ans = 3 >> rank([a b]) %计算增广矩阵[A b]的秩 ans = 3 非齐次线性方程组系数矩阵A的秩为3，增广矩阵的秩为3，等 于未知数个数3，方程组有唯一解。
取值
最近插值 线性插值，为interp1的缺省方法 分段三次样条插值 分段三次Hermite插值 同pchip
不同方法插值的结果

在[0,2π]区间内生成11个等距的离散点， 计算函数y=sin(x)的数值
最近插值 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0 1 2 3 4 5 6 7
插值的数学定义


已知由g(X) (可能未知或非常复杂)产生的n+1 个离散数据(xi, yi), i=0,1,2,…,n，且这n+1个互 异插值结点满足a=x0<x1<x2<…<xn=b，在插值 区间[a,b]内寻找一个相对简单的函数f(x)，使 其满足插值条件f(xi)=yi, i=0,1,2,…,n。再利用已 求得的 f(x)计算任一非插值结点x*处的近似值 y*=f(x*)。其中f(x)称为插值函数，g(x)称为被 插值函数 从计算的观点看，插值就是用一个简单函数在 某种误差范围内近似的代替原目标函数关系式
18.1741
视频 演示
例7-5计算结果的比较
x 插值方法 2.55 2.63 2.77 2.86
真实值
12.8071
13.8738
解：在Matlab中输入 >> a=[-1 -2 4; 2 1 1; 1 1 -1]; >> rank(a) ans = 2 齐次线性方程组系数矩阵A 的秩为2，小于未知数个数3，方 程组有无穷多解

;不能少
计算系数矩阵A 的秩
例7-2—（2）
7 x1 28x3 1 (2) 28x2 x3 39 28x 196x 7 3 1
2000年，MATLAB 6.0版
2001年，MATLAB 6.1版 2002年，MATLAB 6.5版 2004年，MATLAB 7.0版
Matlab 的优点



语法简单易学，编程效率高 高质量、高可靠的数值计算能力 强大的矩阵运算能力 高级图形和数据可视化处理能力 提供600多个常用算法内建函数，以及众 多面向应用的工具箱
2.3 线性方程组直接求解

例7-3 求以下方程组的解
1 5 x1 6 x2 x 5x 6 x 0 1 2 3 x2 5 x3 6 x4 0 x3 5 x4 6 x5 0 x4 5 x5 1

步骤 a 判断解的情况
1 Matlab 基础知识
1.1 Matlab 简介



1967年由Clere Maler用FORTRAN语言 设计和编写 1984年Mathworks公司用C语言完成了 Matlab的商业化版本并推向市场 经过20余年的改进，Matlab已发展成为 一个具有极高通用性的、带有众多实用 工具的运算平台，成为国际上广泛认可 的优秀科学计算软件

y y1
y=?
y0
y1 y0 y y0 ( x x0 ) x1 x0
x0
x
x1
x
并求取在该点的函数值
线性插值的优点是简单，快捷， 特别是对于插值结点间距较小的 情况可以取得令人满意的精度
3.2.2 二次插值
y


又称拉格朗日三点差值 根据三个已知点 [x0, y0], [x1,y1] , [x2,y2] (x0<x1<x2)， 构造二次多项式插值函数 y=a0+a1x+a2x2，并用该函 数计算在x处的y值 二次插值公式
b1 矩阵除法，验证解 b2 逆矩阵法 b3 rref
例7-4

求下列方程组的解
x1 x2 x3 x4 1 x1 x2 x3 x4 1 2 x 2 x x x 1 2 3 4 1
视频 演示
3 数据插值
3.1 数据插值简介

常用语法：YI = interp1(X,Y,XI,'method')

式中

X,Y───为已知数据点的x, y值； XI─── 为待插值数据点的x值； YI─── 为返回的插值结果； method───用于指定所采用的插值方法
插值函数interp1提供的插值方法
项目
nearest linear spline pchip cubic

在工程领域，许多实验数据常以列表函数或表 格的形式存在，如水黏度随温度的列表函数
t/℃ 0 10 20 30 40 50 60
μ/mPa·s
பைடு நூலகம்
1.788
1.305
1.004
0.8012
0.6532
0.5492
0.4698

在实际使用时，有时需要获得介于表中两个温 度结点之间（如15℃，25℃）的黏度值。而这 些数据未在表中出现，需要我们根据已知的数 据估算出表中未出现的温度点的黏度数值，这 一技术称为插值技术
演示文件 查找帮助 帮助索引 全文帮助
刷新
查找 打印

命令窗口查询


help lookfor
Demos
若不知函数 确切名，可 “Lookfor关 键词”可查

联机演示系统

help
Help全部主题
Help指定函数
例7-1

查找包含“diff”关键词的函数

>> lookfor diff

2.2 线性方程组解的判断

齐次线性方程组AX=0，其解的情况可以 通过系数矩阵A的秩和未知数个数n的关 系来判断



如果系数矩阵的秩为n，方程组只有零解， x=0 如果系数矩阵的秩小于n，方程组有无穷多 解 如果系数矩阵的秩大于n，方程组无解
非其次线性方程组解的情况

非齐次线性方程组AX=b，根据系数矩阵 A的秩、增广矩阵B=[A b]的秩和未知数 个数n的关系来判断其解的情况
最近插值
线性插值 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0 1 2 3
线性插值
4
5
6
7
分段三次样条插值 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0 1 2 3 4
分段三次 样条插值
分 段 三 次 Hermite插 值 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8
3.2 插值方法

线性插值 二次插值 其他插值方法


最近（nearest）插值法 样条曲线（spline）法 埃尔米特（Hermite）法
3.2.1 线性插值
又称两点插值 已知两个数据点[x0, y0], [x1, y1] (x0<x1)，求对应 于x (x0<x<x1)的y值 解法：由x0,y0,x1,y1构造直 线方程
Matlab二维作图
90 120 0.8 0.6 150 0.4 0.2 180 0 South 210 330 30 West 1 60 North
240 270
300 East
1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1
1 data1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1 2 3 4 5 6
1.2 Matlab 的界面
剪切 打开文件 新建文件 复制 粘贴 撤销 重做 启动simulink 导航器 菜单栏 性能改进工具 当前路径 打开目录 打开上一级目录
目录窗口
工作空间 窗口标签 历史窗口 命令窗口 开始按钮
1.3 Matlab 的帮助功能

联机帮助系统
“Help”下拉菜 单中“Full Product Family Help”命令打开 联机帮助系统
1
2
3
4
5
6
7
8
9
Matlab三维作图
40
1
30
0.5
20
0
10
-0.5
0 1 0.5 0 -0.5 -1 -1 -0.5 0.5 0 1
-1 -1
0 -0.5 0 0.5
1
1
-1
0.5 1 0.8 0 0.6 0.4 0.2 -0.5 2 1 0 -1 -2 -2 -1 -5 -5 1 0 0 0 2 0 5 5
用户输入 的命令
查询结果



SETDIFF Set difference. DIFF Difference and approximate derivative. POLYDER Differentiate polynomial. DDE23 Solve delay differential equations (DDEs) with constant delays. DDESD Solve delay differential equations (DDEs) with general delays. DEVAL Evaluate the solution of a differential equation problem. …



如果系数矩阵A的秩等于增广矩阵B的秩且 等于n ，方程组有唯一解 如果系数矩阵A的秩等于增广矩阵B的秩且 小于n ，方程组有无穷多解 如果系数矩阵A的秩小于增广矩阵B的秩， 方程组无解
例7-2 判断方程解的情况
x1 2 x2 4 x3 0 (1) 2 x1 x2 x3 0 x1 x2 x3 0

写成AX=b的一般形式，其中
a11 a A 21 as1 a12 a1n a22 a2 n , as 2 asn x1 b1 x b X 2, b 2 x n bs