nf'?????i?1?1?Xi?1??Ki????1????1??1????Ki???2 (9)
现在再回到COMSOL Multiphysics求根计算。作为一个练习，我们用一般PDE
模式来计算求解。我们可以只载入rootfinder.mph或rtfindgen.mph并加以修改来计算这个问题，但是熟悉COMSOL Multiphysics功能也是非常重要的目的。
表5 闪蒸实例 Model Navigator 选择1D，COMSOL Multiphysics: PDE modes:PDE, general形式 Draw菜单 Specify objects: Line，Coordinates弹出菜单，x：0 1 名称：interval，完成 Options菜单： 输入表4中数据 Constants X1，0.0079 X2，0.1281 X3，0.0849 X4，0.2690 X5，0.0589 X6，0.1361 X7，0.3151 Options菜单： 为(8)式中右侧各项定义表达式 Scalar Expressions t1-X1/(1-u*(1-1/K1)) t2-X2/(1-u*(1-1/K2)) t3-X3/(1-u*(1-1/K3)) t4-X4/(1-u*(1-1/K4)) t5-X5/(1-u*(1-1/K5)) t6-X6/(1-u*(1-1/K6)) t7-X7/(1-u*(1-1/K7)) Physics菜单： 选择域：1和2（按住Ctrl键） Boundary settings 选择Neumann边界条件 采用默认设置q=0, g=0完成 Physics菜单： 选择域：1 Subdomain settings 设定F=1+t1+t2+…+t7; da=0 Mesh菜单： Mesh parameters Solve菜单： Solver parameters Post-processing： Point Evaluation 选择Init选项卡，设定u(t0)=0.5 设定最大基元尺寸为1 点击Remesh，完成 选择Stationary nonlinear求解器求解，完成 选择边界：1，表达式：u，完成 Report窗口，值：0.458509，表达式：u
启动COMSOL Multiphysics，打开模型导航栏。如果你已经打开COMSOL Multiphysics对话框，那么将你的工作空间保存为mph文件或者将命令保存为m文件，然后打开“File”菜单选择“New”。按照表5的步骤建立闪蒸模型。注意本例中有两个新增内容——“Options：Constants”和“Options：Expressions”。在这里可以定义常数，然后在COMSOL Multiphysics中任何可以使用纯数字的输入框中使用定义的常数。表达式和常数类似，可以在任何COMSOL Multiphysics数字输入框中使用，但是不同之处在于表达式依赖于因变量。定义的常数和表达式同样可以在后处理过程中使用。后处理过程显示t1和t2为：
Value: -0.013865, Expression: t1, Boundary:1 Value: -0.203441, Expression: t2, Boundary:1
另外求解器日志中经常包含有用的信息。打开“solver”菜单，在最底部选择“View Log”，会弹出求解器日志对话框。它显示了求解器何时开始计算，执行了什么命令，如何从初始状态得到计算结果。这里经过三次迭代，绝对误差达到10－9。如果你的解不是不收敛或收敛很慢的话，很容易得到这样的信息。 练习
1.1 计算方程f(u)?u3?3u2?52u?12?0的根。由于该函数是三次多项式，所以存在解析的无理数解，
u?1u?1?12u?1?12。
u exp(u)-1
1.2 计算方程f(u)?ueu?1?0的根。这是个超越方程，意味着不存在解析的无理数解。如果你使用系数模式，设c=1来帮助收敛。
3. 方法2：步进法数值积分
数值积分是数值分析的支柱。在有计算机之前，科学计算的首要工作是制作特殊方程手册，其中几乎包括了所有特殊常微分方程的解。那么用什么方法计算得到的？就是用一维数值积分。
一维数值积分分为两种：初值问题（IVP）和边界值问题（BVP）。我们下一节再介绍边界值问题。最容易积分的是初值问题，因为所有的初始条件都在一点定义，可以根据一阶导数直接步进积分。显然，如果常微分方程是一阶的，例如：
dydt?f(t)yt??t?yt??tf(t) (10)
那么当△t→0时，(10)中第二项严格成立。它满足欧拉法的条件，是积分一阶常微分方程最简单的办法。对于一维情况，可以简单地根据f在(xn,yn)点的导数向前步进积分，这里n是指积分过程中的第n次离散化步骤。所以，
yn?1?yn?hf(xn,yn)xn?1?xn?h (11)
这里假设导数变化不超过步进量h，这只对线性函数才严格成立。对于任何有曲率的函数，该假设都不成立。下面考虑在图1中，如果步进量h取得过大会造成什么错误。显然，欧拉法需要改进的一点就是增大步进量，因为它需要小步进量以保证准确性。欧拉法也被称作“一阶”准确性，因为误差只能降低到一阶h的程度。
图1 数值积分中的欧拉步进
龙格库塔方法
所以如果我们希望增大步进量，那么就需要一个“更高阶的方法”，随着步进量的减小，误差快速降低。k阶方法误差大约在hk量级。假设我们忽略了曲率，可以通过计算斜率f(x)在xn和xn+1之间几个中间点的值来计算y(x)的曲率。可以首先根据初始导数计算间隔中点的值，然后再用整个间隔中点的微分值来获得二阶准确性。
k1?hf(xn,yn)11??k2?hf?xn?h,yn?k1?22??yn?1?yn?k2?O(h)3 (12)
这样我们通过对两个函数的计算又得到一阶准确性。所以，如果使用一阶方法，
N次计算得到O(1/N)误差，而二阶方法2N次计算得到O(1/4N2)误差，如果使用一阶方法的话，这需要N2次计算。 更高阶龙格库塔方法
我们可以做的更好吗？显然，如果用三个中间点可以获得三阶准确性，用四个中间点可以获得四阶准确性等等。不过更高阶的方法需要更多的编程工作，所以也要考虑我们的时间。但是从本质上讲，我们计算函数的k阶导数不一定光滑。可能在增加“近似准确性”的同时，高阶导数的整体误差会变得很大，如果是这样，每次连续步进都可能导致误差的急剧增大。这表明高阶方法没有低阶方法稳定。
常微分方程通常选用四阶龙格库塔方法积分，这样程序比较紧凑，准确性高，也有很好的稳定性。 其它方法
在数值积分编程的时候，还需要注意两个问题：
数值稳定性：即使使用高阶准确性的方法，你的积分结果与检验值偏差仍然很大，这可能就是因为你的数值方法不稳定。你可以通过降低步进量来获得数值稳定性，但是这也意味着更长的计算时间。如果你的计算量很大，计算时间太长，可以选用半隐式方法，例如预测校正法等。
刚性系统：物理机理起作用时，刚性系统通常会有两个完全不同的长度或时
间尺度。刚性系统可能会出现上面提到的“系统非稳定性”现象，或者非物理现象的振荡。可以参考Gear[2]的书来解决这个问题。
3.1 数值积分简单实例
高阶常微分方程通常通过降阶和步进法求解，考虑如下方程：
dydx22?q(x)dydx?r(x) (13)
除非q(x)和r(x)是常数，那么你很幸运能够从大多数分析方法教材中找到答案。也有一些特殊的q(x)和r(x)可以求得解析解，但是最好能计算出各种情况下的数值解。为什么呢？因为为了搞清楚y(x)，你需要画出曲线，得到解析解时你也需要花费一定的精力在作图上。那么应该怎么做？首先将上面的二阶系统降低为两个一阶系统：
dydxdzdx?z(x)
?r(x)?q(x)z(x)上式的常微分方程能够仿照(11)和(12)进行时间步进法数值积分。举个简单例子：
dudt22?u?0 (14)
降阶产生了两个一阶常微分方程：
du1dtdu2dt??u2
?u1 (15)
假设初始条件为u1=1，u2=0，可以建立零维空间系统来积分这一对常微分方程。
打开COMSOL Multiphysics，等待模型导航栏窗口。如果已经打开一个COMSOL Multiphysics窗口，把工作空间保存为mph文件或者把命令保存为m文件，然后在“File”菜单中选择“New”选项。按照表6的步骤进行设置。
表6 时间步进积分的简单实例 选择1D，COMSOL Multiphysics: PDE modes: PDE, generalModel Navigator 模式，插入因变量：u1 u2 Specify objects: Line。Coordinates弹出菜单，x：0 1 Draw菜单 名称：interval，完成 选择域：1和2（按住Ctrl键） Physics菜单： 选择Neumann边界条件 Boundary settings 保持默认选项q=0, g=0完成 选择域：1 Physics菜单： 设定F1＝－u2，F2＝u1 Subdomain settings 选择Init选项卡，设定u1(tw)=1 Mesh菜单： Mesh parameters Solve菜单： Solver parameters Post-processing： Cross-section plot parameters 这个例子给了我们两个独立变量u1、u2和空间坐标x。注意“Subdomain settings”窗口左上角方程中的矢量标志。由于是矢量变量，所有输入选项卡都是矢量(F)或者矩阵(da)输入。
MATLAB命令linspace(0,2*pi,50)生成50个基元的矢量，给出从0到2π的数据，u1(t=2π)=1.004414。假设解析解是u1(t=2π)=1，结果不够准确（0.4％误差）。之前FEMLAB允许用户在MATLAB和FEMLAB中选择针对时间积分的预置的求解器。也允许用户在求解器变量对话框的通用选项卡中调整容错值（相对和绝对）。将相对误差调整为0.001，绝对误差调整为0.0001，得到一个相对较好的计算结果u1(t=2π)=0.998027。注意累积的全局误差与求解方法准确性0.001的量级一致。
设置最大基元尺寸为1 点击Remesh，完成 选择time dependent求解器，在General选项卡中输入 Times: linspace(0,2*pi,50) 求解，完成 Point选项卡，接受u1的默认设置 General选项卡，接受所有时间的默认设置 完成
图2 一个周期的u1(t)