微型机械液晶引流计算机辅助计算软件设计本科毕业设计论文.docxVIP

机电工程学院毕业设计说明书设计题目: 微型机械液晶引流计算机辅助计算 软件v1.1 设计 目录1 引言11.1 课题简介11.2 课题研究的国内外现状21.3 软件设计方案32 课题相关技术简介52.1 vb简介52.2 电场驱动下液晶引流相关微分方程92.3 数值分析之常微分方程解法102.3.1 euler法102.3.2 龙格-库塔法122.3.3 中心差分法(central difference method)原理133 软件设计流程173.1微分方程分析173.2 龙格-库塔算法在此微分方程中的应用理解183.3 程序设计183.3.1变量定义183.3.2 核心程序段经典四阶龙格-库塔法说明204 设计成果及展望304.1 数据分析304.2 程序设计展望354.3 设计总结41致谢43参考资料441 引言1.1 课题简介本 课 题 来 自 于 机 电 工 程 学 院 刘 春 波 老 师 的 国 家 自 然 科 学 基 金 项 目 即 关 于 建 立 在 液 晶 引 流 基 础 之 上 的 微 流 体 全 新 驱 动 方 式 以 及 控 制 原 理 等 方 面 的 研 究 。 液 晶 引 流 效 应 是 指 在 外 加 磁 场 或 电 场 作 用 下 ， 液 晶 流 动 与 其 内 部 分 子 指 向 向 量 的 排 列 之 间 相 互 作 用 所 产 生 的 一 种 现 象 ， 其 中 内 部 分 子 指 向 矢 指 液 晶 内 部 一 点 附 近 的 分 子 平 均 指 向19。 电 场 驱 动 液 晶 引 流 的 原 理 即 是 在 电 磁 场 的 直 接 作 用 下 ， 指 向 列 液 晶 由 于 引 流 效 应 的 影 响 就 会 产 生 流 动 ， 液 晶 的 表 观 粘 滞 系 数 （ 表 观 黏 滞 系 数 的 计 算 方 法 是 ： 单 位 面 积 所 受 切 向 力 除 以 速 度 梯 度 ） 受 到 这 种 方 式 的 流 动 影 响 ， 从 而 破 坏 液 晶 的 显 示 效 果 , 向 列 相 的 液 晶 和 一 般 液 体 的 主 要 特 征 区 别 在 于 其 分 子 形 状 为 长 棒 形 ， 而 且 它 的 长 轴 所 指 的 方 向 总 体 上 大 致 相 同 ， 一 般 我 们 将 这 一 方 向 定 义 为 液 晶 指 向 矢 的 方 向 ， 在 外 加 电 磁 场 的 作 用 下 ， 指 向 矢 就 会 发 生 连 续 的 变 化 ， 从 而 导 致 液 晶 流 动 的 发 生 ， 因 此 当 在 研 究 液 晶 流 动 的 过 程 中 ， 我 们 需 要 同 时 解 决 两 个 场 ， 一 个 是 分 子 的 配 向 场 ， 也 就 是 液 晶 指 向 矢 所 指 方 向 的 变 化 场 ， 另 外 一 个 就 是 液 晶 流 动 的 流 场 2 9 。 这 两 个 场 相 互 作 用 ， 相 互 影 响 ， 是 我 们 研 究 的 关 键 所 在 。建 立 在 液 晶 引 流 基 础 之 上 的 全 新 微 流 体 驱 动 与 控 制 方 式 是 将 液 晶 弹 性 形 变 理 论 和 l e s l i e - e r i c k s e n 理 论 相 结 合 ， 对 两 平 行 板 内 电 磁 场 引 起 的 液 晶 流 动 性 进 行 研 究 1 9 。 各 变 量 之 间 的 相 互 关 系 在 液 晶 流 动 的 过 程 中 相 当 复 杂 ， 当 我 们 运 用 数 值 分 析 的 计 算 方 法 对 其 进 行 分 析 时 ， 其 运 算 量 相 当 之 大 ， 因 此 我 们 必 须 借 助 计 算 机 等 辅 助 工 具 进 行 计 算 。 本 课 题 是 是 根 据 电 场 驱 动 液 晶 引 流 的 一 些 原 理 ， 来 分 析 电 磁 场 作 用 下 条 状 液 晶 的 一 些 相 关 性 质 ， 而 本 毕 业 设 计 所 研 究 的 微 分 方 程 正 是 在 此 定 性 分 析 中 推 导 出 来 的 ， 然 后 根 据 微 分 方 程 的 具 体 形 式 和 它 的 初 始 计 算 条 件 ， 选 择 合 适 的 数 值 计 算 方 法 进 行 计 算 ， 并 输 出 保 存 计 算 数 据 ， 为 绘 图 做 数 据 储 备 。刘老师在项目分析中，已经推导出了电场驱动液晶引流的相关微分方程，并且导师推荐了进行计算机仿真计算的方法即四阶龙格-库塔（runge-kutta）算法。本课题最终的目的是为了开发一个计算机仿真计算软件，对电场驱动液晶引流的相关计算公式（多元微分方程组），利用数值分析中的龙格-库塔算法和中心差分法进行计算机仿真计算，在我们计算前已经设定好了一定的初始计算条件、常数项系数以及时间步长和空间距离步长，我们用四阶龙格-库塔算法进行迭代计算，可以计算出一定时间段内每一时刻的近似函数推算数据，然后将计算结果以txt文件形式进行保存，对数据计算中出现的的误差和问题进行分析，也就是要进行数据检测。空间上可以采用中心差分法，时间上可以采用龙格-库塔法进行计算，要求计算精度越高越好，至少能达到。1.2 课题研究的国内外现状液 晶 分 子 在 显 微 镜 下 呈 长 条 状 ， 在 电 磁 场 的 作 用 下 ， 液 晶 由 于 引 流 效 应 的 影 响 会 产 生 流 动 ， 近 年 来 微 流 体 驱 动 技 术 已 经 取 得 了 飞 速 的 发 展 ， 研 究 员 们 越 来 越 留 心 注 意 微 米 乃 至 纳 米 尺 度 构 件 中 流 体 的 驱 动 与 控 制 技 术 ， 同 时 液 晶 独 特 的 流 动 特 性 也 成 为 了 研 究 的 热 点 ， 诸 多 方 向 的 研 究 也 方 兴 未 艾 ， 最 初 这 些 研 究 主 要 集 中 于 液 晶 流 动 的 速 度 场 与 配 向 场 相 互 作 用 及 其 e r m r 特 性 研 究 等 方 向 ， 但 是 这 一 阶 段 的 研 究 局 限 于 从 流 体 力 学 的 角 度 探 索 液 晶 流 动 的 性 质 ， 还 是 把 液 晶 作 为 普 通 流 体 的 一 种 1 9 。 2 1 世 纪 以 来 ， 纯 理 论 的 研 究 开 始 向 应 用 方 向 转 变 ， 首 先 出 现 了 利 用 引 流 效 应 驱 动 液 晶 内 微 纳 米 级 粒 子 的 研 究 ， 但 由 于 粒 子 的 存 在 影 响 了 液 晶 的 流 动 特 性 ， 而 且 驱 动 电 压 很 高 ( 1 0 0 v ) ， 这 就 使 理 论 计 算 与 实 验 研 究 的 结 果 由 于 自 由 电 子 的 注 入 相 差 较 大 。 从 液 晶 自 身 特 点 出 发 ， 要 实 现 精 确 驱 动 ， 必 须 考 虑 液 晶 分 子 特 殊 形 状 及 动 作 特 征 的 液 晶 流 动 机 制 及 影 响 因 素 。 在 静 态 理 论 中 ， 最 精 确 的 连 续 体 模 型 是 由 o s e e n 与 z c h e r 提 出 的 ， 该 模 型 充 分 考 虑 了 向 列 相 液 晶 分 子 独 特 的 长 棒 状 特 点 ， 后 来 经 过 f r a n k 的 补 充 形 成 了 目 前 为 止 一 直 沿 用 的 曲 率 弹 性 形 变 理 论 ， 学 术 界 最 具 影 响 力 的 动 态 理 论 是 由 l e s l i e 和 e r i c k s e n 联 合 提 出 的 l e s l i e - e r i c k s e n 理 论 （ l - e 理 论 ） ， 分 别 从 动 量 、 质 量 、 能 量 和 角 动 量 的 角 度 阐 明 了 外 力 对 于 向 列 相 液 晶 配 向 场 和 流 场 的 影 响 ， 并 给 出 了 适 用 于 向 列 相 液 晶 的 本 构 方 程 ， 以 这 两 种 理 论 为 基 础 ， 对 两 平 行 板 内 施 加 电 场 后 所 引 起 的 液 晶 流 动 正 在 研 究 当 中 ， 其 主 要 是 为 了 揭 示 流 动 形 成 的 机 制 及 影 响 因 素 1 9 。 目 前 关 于 微 流 体 和 液 晶 的 研 究 ， 在 日 本 ， 较 之 于 世 界 其 他 各 国 ， 成 就 斐 然 ， 处 于 领 先 地 位 。目前关于微流体和液晶的研究，在日本，较之于世界其他各国，成就斐然，处于领先地位。1.3 软件设计方案 本方案是在windows xp系统visual basic6.0环境下编译的。根据visual basic语言的特性,具体的设计步骤如下： 1）设定a开头的变量代表变化量，b开头的变化量代表已知量，c开头的变化量代表暂确定量。 2）定义d开头的代码变量表示在程序运行中使用到的过程变量。 3）定义保存数据所用变量。 4）把计算初始已知的数据和条件赋给设定好的常系数变量。 5）设定相关数组的区间大小。 6）给程序运行的循环变量设定最初数据。 7）设置保存文本的具体格式。 8) 核心计算代码段，即四阶runge-kutta法，分五步完成： 1. 计算k1 2. 计算k2 3. 计算k3 4. 计算k4 5. 利用runge-kutta合成计算说明：计算过程由于用到龙格-库塔算法和中心差分算法，根据二者算法的特性，又由于微分方程是含有两个未知数、五个因变量的方程组，计算前先将定义区域划分网格，网格的大小即为时间和空间距离的步长，计算每一时刻在每一空间位置处相应的函数值，使得生成的数据可以在绘图软件中绘制成动态变化图像。k1、k2、k3、k4的计算有着相似的计算步骤，但是应注意关于因变量在循环计算过程中的中间迭代。 9）保存这一时刻对应空间距离的函数值。对应划分好的的每一时刻都有对应的一组函数值，在坐标系中可以绘出相应的图像。 10）预估程序运行总耗时量。 2 课题相关技术简介2.1 vb简介 basic是（beginners all-purpose system instruction code）的缩写，具有面向普通使用者和易学易用的特点，于1964年创建，在20世纪70年代末得到了很大的发展，它作为一种古老的程序设计语言对计算机的普及和推广起到了举足轻重的作用。basic语言几乎是拥有用户最多的计算机语言。visual的英文翻译是“视觉的”或是“可看得见的”，在此指开发图形用户界面(graphical user interface,gui)的方法，即可视化程序设计。这种编程方法在描述界面外观的位置时不需要大量代码，而只要把预先建立的控件，像使用“画图”之类的绘图程序那样绘到屏幕上即可13。 visual basic因其易学易用，受到广大程序开发人员的喜爱。visual basic的编程语法结构与qbasic基非常相似，凡是学过qbasic的人很容易掌握visual basic。我们在大学阶段通常学习的是c语言或是c+语言，他们属于高级的编程语言，相对抽象，然而vb却更接近正常语句，相对简单，因此研究员们在开发一些应用软件时会首先选择vb语言。 visual basic是当今全球最畅销的编程语言之一，它简单易用，是一种可视化的语言，被众多软件开发者所青睐，有以下特点：1）.可视化的编程 在以前我们编写c 程序或者是c+程序时，我们必须且只能是在编写代码过程中，将软件相应的界面设计并体现在呈现在代码过程中，如果软件界面设计不妥，只能是在程序代码运行调试之后方才能够发觉，如果要修改只能关闭软件，然后返回程序代码段，修改代码，再次调试运行，反复上述操作直至软件界面达到用户要求为止，如此显然增加了诸多不变，柔性降低，不利于实现软件设计的灵活机动性。而vb是将相应的程序段代码比如文本输入，设定并整合为一体，形成控件，以图标的形式排列在工具栏处，编程前用户可以先将控件排列布局，形成美观合乎要求的界面后，然后在进行相应的程序编制，简单而灵便。2）.面向对象的程序设计 visual basic是面向对象程序设计语言。面向对象的程序设计方法，是指把程序和数据封装作为一个实体，程序设计针对这些对象进行，不必重复编写大量的代码。3）.机构化程序设计语言 visual basic是属于相对比较高级的编程语言，并且代码相对通俗，更接近人们说话的方式，对于非编程专业人员理解起来相对容易，在vb中生成的软件，在脱离原编程电脑后，对于另外一台电脑只要是windows环境下，基本上都可以运行，这是vb的优越之处。4）.事件驱动编程机制 visual basi程序设计中对对象的操作要通过事件来完成，一个对象可对应多个事件，一个事件要通过一段程序来执行。例如visual basic窗体对象中有load事件，对应一段程序实现指定的操作。5）.访问数据库visual basic系统有很强的数据库管理能力，利用数据控件和数据库管理窗口，可以直接之间建立或处理access格式的数据库；同时，visual basic还能编辑和访问外部数据库，如foxpro、dbase等；visual basic还提供开放式数据链接（odbc）功能，通过它可以访问和连接后台大型数据库，如sql server、oracle等。一般情况下，我们操作vb的流程如下：1. visual basic的启动 1）在桌面上依次选择”开始”“程序”“microsoft visual basic中文版”程序组“microsoft visual basic中文版”，从而启动vb，屏幕上显示“新建工程”对话框，如图2-1所示 图2-1 启动后“新建工程”对话框 2）双击上述对话框中“标准exe”或直接单击“打开”按钮，进入vb开发环境，如图2-2所示 图2-2 visual basic 集成开发环境2. 编写相关的程序代码设计界面时，根据自己的需要，在图2-2中form1窗体中添加相应控件并布局（在vb中空间是软件中封装起来的应用程序，是可视化特征存在的根本，它可以很方便用户不必要在代码中设置软件界面，可以根据自己的喜好随拉随放，自由排列），然后可以在右侧设置相应控件的基本属性，当双击相应控件就可以进入代码编辑的窗口，然后单击“选择对象”下拉列表框右侧箭头，并从中选择“form”对象，再从“选择事件”下拉列表框中选择“load”事件，此时在代码窗口中将出现事件过程的框架。我们就可以将程序代码输入窗口中。3 保存工程文件点击执行“文件”下拉菜单“工程保存”的命令，或者单击工具栏上的“保存”按钮，则vb系统就会将所有的有关内容按一定路径保存。如果这是系统第一次保存该文件，vb就会出现“文件另存为”的对话框。用户可以根据自己的意愿设置待保存的文件名，并指定相应位置和保存路径。如果这不是系统第一次保存文件，则系统将会直接地将工程中的所有文件以原文件的名称进行保存。如果我们要以新的文件名进行保存，那么可以从文件菜单中单击“工程另存为”命令，同样可以出现文件另存为对话框，用新的文件名保存此工程文件。4. 运行调试程序单击工具栏中“运行”下拉选项中“启动”命令，系统将会进入运行状态，如果所编写程序代码中没有出现任何错误，则该程序将继续执行；如果程序代码中出现错误，程序代码运行出错的消息框将会出现在窗体中。如图2-3所示 图2-3 程序运行出错提示这时用户点击确定，然后单击“运行”下拉选项中“结束”按钮，则会结束正在运行的程序，并回到编写程序的模式，用户就可以从代码窗口中修改运出现错误的程序代码。或者是单击“调试”下拉选项中“逐语句”，进入程序调试状态，此时界面中将会出现程序代码窗口，运行中有错误的程序代码行上将会以黄色作为提示，黄色光标将会停留在出错的程序段上 。修改错误完成后后，就可以继续按原来的步骤运行。如此反复直到没有错误为止。2.2 电场驱动下液晶引流相关微分方程 在课题研究中刘春波老师已经推到出了电场驱动下液晶引流的相关微分方程，具体如下：其中已知量如表所示(pas) (n) (f/m ) 1020 -8.6 -0.4 8.9 5.9 -3.110126.37 3.81 8.6015.7 5.7 1000 初始值：u=0 ; w=0; nx=cos1; ny=cos1; nz=0; e=5;条件设定：2.3 数值分析之常微分方程解法数值分析是一门研究各类数学算法的学科，可以帮助寻求我们在进行科学研究中所遇到的数学问题近似解，他的基本思想应属于高等数学中的极限理念，也就是满足一定的精度条件下，无限逼近所求的数学问题的真值，这就好比机械加工中的公差概念，只要产品在一定的公差范围内，结果都是满足要求的，他们的基本思想是相通的。常见的算法有求解非线性方程的二分法、牛顿发，多项式插值中的三次样条插值法、hermite插值法，求解常微分方程中的euler法、龙格-库塔法，在此我们主要谈一下求解常微分方程中的euler法、龙格-库塔法。2.3.1 euler法对于给定初始值的一阶导数的微分方程 y=f(x,y) （cxb） y（a)=c 我们通常假设式在区间c , b上存在唯一解y(x),且y(x)在c , b上具有所需的光滑性，也就是y(x)在c , b上存在所需阶数的连续导数。为求式的数值解，我们采用离散化方法在区间a , b上取一组节点： c=b称h=（0in-1）为步长，简便起见常规定其为定值 =a+i*h （0in，h=（b-c）/n）离散化方法的基本特点是按某一递推公式，按节点从左至右的顺序依次求出y（）的近似解取=c，如果计算只用到之前一步的，则称该方法为单步法，如果用到前r步的值，就称作r步法。euler法是解决这类初值问题的常用数值方法，将式的两端在区间，上积分，就得到dx=f（x,y(x)）dx 即 y（）=y()+f（x,y(x)）dx 由泰勒公式并同时消去无穷小量，可得到： =+h*f（,）由于=c，这样就能够依次求出，这种方法就称作euler法。 euler法的几何意义是对于比如二维平面下的一条曲线，当按照欧拉算法，就可以求出一系列坐标点x（0）、 x（1）x（n-1）、x（n）所对应的近似解、，相邻两点彼此连接，就可以逼近拟合原曲线，当步长足够小，算法恰当就可以达到一定的精度，所以euler法也称折线法。对于 式，应用积分中值定理可得 y（）=y()+ h f（+h,y(+h)）此处f（+h,y(+h)）称作，上的平均斜率，记作，即 =f（+ah,y(+ah)）可以看出，只要为平均斜率提供一种算法，就可以得到一个微分方程的数值计算公式，上述euler公式仅取一个点的斜率值f（,）作为平均斜率的近似值，因此精度相对较低。然而改进的euler公式利用了两个点的斜率值=f（,）与=f（,）的平均值作为平均斜率的近似值，即 1/2*（+）因此改进欧拉公式为 ： =+h*1/2（+） =f（,） =f（,） =c 2.3.2 龙格-库塔法 在上述分析中可知，由泰勒公式或者积分特性已经可以对微分方程进行近似求解，即将微分方程的形式转化成代数方程组，但是我们在长期使用过程中发现，其计算精度相对较低，面对复杂而工作量庞大的计算任务时，无法达到预期的精度，即便用改进的euler算法，精度依然受限。因此数学家们在长期的研究中总结出了runge-kutta算法，其推导思想近似于euler算法，当是单步计算时，也就是一阶runge-kutta算法和euler算法相同，即 (此处h称作步长，且h=xi+1-xi，即相邻离散点的距离，为方便计算，一般令其为常值，整体在定义区间的含义是足够小等间距分割。)由上述 euler算法介绍可知，当在相邻两点之间插入一点进行斜率计算，并对两斜率加权平均后可以使相应的计算精度大大提高，由此我们同样得出了二阶runge-kutta算法的计算公式，只是推导过程不同，注意区别，具体相见文献30。因此二阶runge-kutta算法中点公式为： =+h* 同样的道理，我们依此类推，经数学推导、求解，可以逐个增加斜率，由于三阶runge-kutta算法在工程计算中不常使用，我们不再做介绍，当在相邻两点之间插入四点，求解四个关联的斜率时，同理可以得出四阶runge-kutta算法公式，也就是在数学以及其他领域经常使用的经典四阶runge-kutta算法： 一般情况下，随着所插入的点数的增加，也就是所求斜率数量的增加，近似计算的函数值精度将会随之增加，但是四阶龙个-库塔以下的算法精度达不到四阶，随着计算总量的增加，误差的逐渐累积，会使后续数据失真，甚至出现紊乱，丧失参考价值。2.3.3 中心差分法(central difference method)原理差分的基本概念：设x的解析函数y=f（x）,函数y对x的导数为 这里dy、dx分别为因变量和自变量的微分形式，dy/dx为函数对自变量的导数，又叫微商，在上式中为函数对自变量的差商，、分别为因变量和自变量的差分。中心差分的方程为: 中心差商的方程为: 中心差分法属于有限差分法的一种，在求解流动计算方程时（此时方程中有两个自变量，定义区间为二维平面，属于求偏微分）的一般步骤是： 图2-4 时间和空间的网格划分1. 先将求解区域划分为差分网络，即利用有限个离散点代替连续的求解区域，如图2-4所示。2. 将要求解的物理量，例如：密度、速度等分配在相应位置的网格点上。3. 将偏微分方程的各微分项用差商代替，从而也就将整个微分方程组转化为代数方程组的差分方程，这样也就得到了各离散点的差商方程组。4. 迭代计算出差商方程组的所有解，即各网格点处未知待求变量的数值解。 由上述关于差分和差商定义的分析可知，当自变量的差分，即增量，无限趋近于0时，就可以由差商得到导数，因此在数值分析计算中，常用差商代替求解导数值。中心差分法是建立在有限差分的基础之上，也就是利用数值的方法分割离散点求解关于一阶和二阶导数的问题，这里用位移对时间的一阶和二阶导数来加以说明，其一阶导数为速度，二阶导数为加速度。其基本思想是把速度和加速度向量用只含有位移的某种方程组合来表示，舍去偏微分、求导等复杂的数学形式，也就是把微分的计算形式转化为代数方程组的形式进行数值计算，并在一定时间区间内求得每个短时间段内下一时刻对上一时刻的的迭代计算公式，进而就可以计算出整个时间区间内各相应位置处的速度和加速度。推导如下：由泰勒公式： (xx+h) 即 又 (x-hx) 即 +且当h0时，二阶导数项为高阶无穷小，舍去可得： 同理仍用泰勒公式推导计算x处二阶导数可得： 以位移-时间图像为例，如图2-5所示来说明用中心差分法计算一阶导数即速度，用u代表位移，t1时刻位置为u(t1)，t2时刻为u(t2)，t3时刻为u(t3),由中心差分法可知，当t3-t2=t2-t1=，且当为极小值时，在一定精度内，速度的数值计算方法为 以此方法进行递推，用u(i-1)表示物体在时刻的位置，u(i)表示物体在时刻的位置，u(i+1)表示物体在时刻的位置，其中-=-=，t为时间步长常数，如此速度v的中心差分为： 加速度a的中心差分为：速度和加速度的微分形式分别为 v = u(i) = a = u(i) = 很显然我们把微分形式转化成了代数方程的形式，对于我们实现计算机vb语言编程，提供了可能。 图2-5 位移-时间坐标系 3 软件设计流程3.1微分方程分析 在此五个微分方程中除去常数系数之外，总共含有五个因变量，即u、w、nx、ny、nz,两个自变量，即时间t，空间距离y，u、w、nx、ny、nz为t、y和函数值构成的三维空间的变量，其中nx、ny、nz要对时间t进行一次求偏导，对y进行一次和二次求偏导，(见上述五个微分方程)，针对每一个方程可以理解为：在任一时刻t，其任一因变量相对于自变量即空间距离y，近似为二维曲线，亦即把时间变量看作常量，每个方程即为因变量对y的二维曲线。(此为程序设计核心部分，即总循环结构的逻辑)编程初步难点：1. 因变量对自变量y的一阶和二阶导数以及t的一阶导数如何同时体现出来，即和等对策：鉴于对每一时刻t，因变量对自变量是二维曲线，可以先将时间离散，对每一时刻t固定，然后逐步依次对y求函数值，待求至y取值范围终点，在返回求下一时刻t依次关于y的函数值，如此延续，并将计算数据输出。1. 形如()如何使用中心差分法表示。对策：分别表示，然后合成，对于的表示比较容易，即对任一空间位置=(u(i+1)-u(i)/y，那么()可以把作为一个整体进行中心差商，由于nx、ny的初值是已知的，那么就可以依次在各离散点求其中心差分值，如此二者合成即可。3.2 龙格-库塔算法在此微分方程中的应用理解 由于空间上使用中心差分法，时间上使用龙格-库塔法，对应于公式其f函数相当于因变量u、w、nx、ny、nz对时间t的一阶导数，即、，而四阶龙格-库塔算法中k1、k2、k3、k4中yi、yi+1/2*h*k1、yi+1/2*h*k2、yi+h*k3均等价于一阶导数的自变量，在程序设计中应加以注意，例如 对于 y=-2x (0x1.2) y(0)= 1 利用四阶龙格-库塔算法公式得： k1=-2* k2=-2*(+h/2)(+h/2*k1) k3=-2*(+h/2)(+h/2*k2) k4=-2* 计算出来后对四个斜率值加权平均带入迭代公式 即可逐步求得结果。 综上所述，数值分析的基本思想是对于连续的函数曲线，采用有限的离散点分割，然后利用一个算法求取某一点的近似函数值，当步长足够小，算法恰当，可以达到相当高的拟合精度，简而言之即用离散点去逼近拟合函数曲线。3.3 程序设计3.3.1变量定义 a开头代表变化量，b代表已知量，c代表暂确定量，d代表中间变量活循环变量。说明： a开头变量主要是针对时间t的迭代循环，形如a_u()是代表二维数组，由于每一时刻的函数是一条曲线，那么就假定a_u(0，i)代表上一时刻即当下时刻的已知函数值，而用a_u(1,i)代表下一待求时刻的函数值，当0时刻的函数值求完之后，求出1时刻的函数值，最后再将1时刻的函数值重新赋给0时刻作为已知值，从而进行下一时刻的运算，如此循环至终。即程序 for j = 1 to n - 1 a_u(0, j) = a_u(1, j) a_w(0, j) = a_w(1, j) a_nx(0, j) = a_nx(1, j) a_ny(0, j) = a_ny(1, j) a_nz(0, j) = a_nz(1, j) next j b开头的变量为程序初始确定已知量，如b_ea 代表微分方程中, b_r1 代表方程中r1, b_r2 代表方程中的r2，由于方程中e是暂定值，因此用c_e代表方程中e，根据情况随时可以改变确定e的值由于保存数据所用，程序在中间过度计算时需要许多中间变量，具体说明：以u为例 (其余四个因变量雷同考虑)，在四阶龙格-库塔算法中， 此处，即为f函数，对于下一个点斜率即k2中+h/2*k1，用d_u_k1表示，依次中+h/2*k2用d_u_k2表示,同理中+h*k3用d_u_k3表示，从而方便进行下一点的斜率计算。相应代码为 d_nx_k1(j) = a_nx(0, j) + dnx_dt_k1(j) * b_dt / 2# d_ny_k1(j) = a_ny(0, j) + dny_dt_k1(j) * b_dt / 2# d_nz_k1(j) = a_nz(0, j) + dnz_dt_k1(j) * b_dt / 2# d_u_k1(j) = a_u(0, j) + du_dt_k1(j) * b_dt / 2# d_w_k1(j) = a_w(0, j) + dw_dt_k1(j) * b_dt / 2#简单理解，此处类似于 f（x）=那么 只不过用z=x+1 代替。3.3.2 核心程序段经典四阶龙格-库塔法说明1. 总循环逻辑 do = calculate_step = calculate_step + 1 t = b_t_start + save_times * save_dt + calculate_step * b_dt if calculate_step = save_step_t then calculate_step = 0 save_times = save_times + 1 loop while t b_t_end =说明：由之前论述，由于每一时刻五个因变量相对于自变量为一条二维曲线，因此程序需要计算出在设定好的各时刻时每一离散点处的函数值，并加以保存。故而计算前需要将总时程进行离散，分割成各个时刻。此处不是靠时间的直接递增来进行下一刻的计算，而是靠计算步次的逐步加一，即calculate_step = calculate_step + 1，计算步次的起始位置是从上一保存时刻之后算起，直到时间递增到设定下一保存的时刻为止，进行下一时刻计算前务必需要将calculate_step 清零，即程序段：if calculate_step = save_step_t then calculate_step = 0因为在各保存时间间隔处数据需要保存，这里定义一个保存次数变量 save_times，每一时刻计算完成自动加1，如是时间的计算方法为 t=开始时间+保存次数*保存时间间隔+计算步次*时间步长即 t = b_t_start + save_times * save_dt + calculate_step * b_dt这样便能够实现时间的逐步递增。1. 计算k1 for j = 1 to n - 1 du_dy(j) = (a_u(0, j + 1) - a_u(0, j - 1) / 2# / b_dy dw_dy(j) = (a_w(0, j + 1) - a_w(0, j - 1) / 2# / b_dy dnx_dy(j) = (a_nx(0, j + 1) - a_nx(0, j - 1) / 2# / b_dy dny_dy(j) = (a_ny(0, j + 1) - a_ny(0, j - 1) / 2# / b_dy dnz_dy(j) = (a_nz(0, j + 1) - a_nz(0, j - 1) / 2# / b_dy d2u_dy2(j) = (a_u(0, j + 1) - 2# * a_u(0, j) + a_u(0, j - 1) / b_dy / b_dy d2w_dy2(j) = (a_w(0, j + 1) - 2# * a_w(0, j) + a_w(0, j - 1) / b_dy / b_dy d2nx_dy2(j) = (a_nx(0, j + 1) - 2# * a_nx(0, j) + a_nx(0, j - 1) / b_dy / b_dy d2ny_dy2(j) = (a_ny(0, j + 1) - 2# * a_ny(0, j) + a_ny(0, j - 1) / b_dy / b_dy d2nz_dy2(j) = (a_nz(0, j + 1) - 2# * a_nz(0, j) + a_nz(0, j - 1) / b_dy / b_dy next j说明：此处是利用中心差分法分别计算0时刻（当下时刻）五个因变量相对于自变量（即空间距离y）的一阶和二阶偏导，为计算k1、k2、k3、k4作数据储备，计算前应先将定义区间（关于时间t和空间距离y）进行网格划分，形如图2-4所示，利用程序循环以代数方程组的形式迭代计算每一时刻的一阶和二阶导数值，此处b_dy代表差分公式中的，a_u(0, j + 1) - a_u(0, j - 1)代表函数u的差分值。 for j = 1 to n - 1 dnx_dt_k1(j) = _ a_ny(0, j) * 2# * (b_k3 - b_k2) * dnx_dy(j) * dny_dy(j) / b_r1 _ - a_ny(0, j) * b_a2 * du_dy(j) / b_r1 _ + (b_k2 + (b_k3 - b_k2) * a_ny(0, j) * a_ny(0, j) * d2nx_dy2(j) / b_r1 式2a（共3项，每项一行） dny_dt_k1(j) = _ a_ny(0, j) * b_ea * c_e * c_e / b_r1 _ - b_a3 * (a_nx(0, j) * du_dy(j) _ + a_nz(0, j) * dw_dy(j) / b_r1 _ - (b_k3 - b_k2) * a_ny(0, j) * (dnx_dy(j) * dnx_dy(j) _ + dny_dy(j) * dny_dy(j) _ - dnz_dy(j) * dnz_dy(j) / b_r1 _ + (b_k1 + (b_k3 - b_k2) * a_ny(0, j) * a_ny(0, j) * d2ny_dy2(j) / b_r1 式2b（共7项，每项一行） dnz_dt_k1(j) = _ a_ny(0, j) * (2# * (b_k3 - b_k2) * dny_dy(j) * dnz_dy(j) _ - b_a2 * dw_dy(j) / b_r1 _ + (b_k2 + (b_k3 - b_k2) * a_ny(0, j) * a_ny(0, j) * d2nz_dy2(j) / b_r1 式2c（共3项，每项一行） next j说明：此处为计算nx、ny、nz的k1值，为后续计算做准备，注意采用的dnx_dt_k1(j)和dny_dt_k1(j)的表达形式，他们都是函数的一阶导数的表达形式同时也是k1值，此处只需将相应的已知变量带入相应已知微分方程即可，后面u、w类似，不再详述，此处留心 a_ny(0, j)，在计算k2、k3、k4时注意函数值的迭代（详见程序代码）for j = 0 to n d_1a_1_ig(j) = d_nx_k1(j) * d_nx_k1(j) * d_ny_k1(j) * d_ny_k1(j) 式1a第1项对y的积分（不含系数） d_1a_2_ig(j) = d_nx_k1(j) * d_ny_k1(j) * d_ny_k1(j) * d_nz_k1(j) 式1a第2项对y的积分（不含系数） d_1a_6_ig(j) = d_ny_k1(j) * dnx_dt_k1(j) 式1a第6项对y的积分（不含系数） d_1a_7_ig(j) = d_nx_k1(j) * dny_dt_k1(j) 式1a第7项对y的积分（不含系数） d_1a_11_ig(j) = d_nx_k1(j) * d_nz_k1(j) 式1a第11项对y的积分（不含系数） d_1b_1_ig(j) = d_nx_k1(j) * d_ny_k1(j) * d_ny_k1(j) * d_nz_k1(j) 式1b第1项对y的积分（不含系数） d_1b_2_ig(j) = d_ny_k1(j) * d_ny_k1(j) * d_nz_k1(j) * d_nz_k1(j) 式1b第2项对y的积分（不含系数） d_1b_6_ig(j) = d_ny_k1(j) * dnz_dt_k1(j) 式1b第6项对y的积分（不含系数） d_1b_7_ig(j) = d_nz_k1(j) * dny_dt_k1(j) 式1b第7项对y的积分（不含系数） d_1b_10_ig(j) = d_nx_k1(j) * d_nz_k1(j) 式1b第10项对y的积分（不含系数） next j 求解各二阶偏导项（不含系数） for j