燃烧矿物燃料对碳循环影响的常微分方程模型

燃烧矿物燃料对碳循环影响的常微分方程模型的求解 燃烧矿物燃料对大气中二氧化碳有何影响 尽管今天二氧化碳 的数量只占大气的大约百万分之 350 但其任何增加都有可能对气 候产生深远的影响 在 LightHouse Foundation 维护的网站 1 上有许 多相关背景的文章介绍 Eric Roden 将 J C G Walker 在 2 中提出的模型引入我们的视线 该模型模拟了许多不同形式碳原子的相互转化 包括存在于大气 浅海及深海的三种状态 模型中五个主要的变量都是时间的函数 p 大气中二氧化碳的分压力 浅海中集中溶解的碳 s 深海中集中的溶解的碳 d 浅海的碱度 s 深海的碱度 d 另外在浅海的平衡方程中还涉及到三个量 浅海中的碳酸氢盐 s h 浅海中的碳酸盐 s c 浅海中气态的二氧化碳的分压力 s p 五个主要变量的变化速率由五个常微分方程给出 大气和浅海之间 的交换涉及到一个固定的特征传输时间 d 和源项 f t 1 tf d pp dt dp s 下面的方程描述了浅海和深海之间的交换 其中 vs 和 vd 分别为两 个区域的体积 21 1 d pp k vdt d s sd s s 1 1sd d d k vdt d 2 1 k vdt d sd s s 1 2sd d d k vdt d 二氧化碳和浅海中的碳酸盐之间的平衡方程由三个非线性代数方程 描述 3 2 1 3 2 2 k k h sssss s 2 ss s ha c s s s c h kp 2 4 模型中涉及到的常量数值如下 64 8 d 2 1 1095 4 2 2 1095 4 12 0 s v 23 1 d v 3 10 4 1 1019 2 k 5 2 1012 6 k 997148 0 3 k 2 4 1079 6 k 源项 f t 描述了现代工业燃烧矿物燃料所造成的影响 时间间隔从一 千年前到未来几千年后 50001000 t 在 t 1000 处的初始值 00 1 p 01 2 s 23 2 d 20 2 s 26 2 d 描述了前工业时代的平衡并保持近似不变 其中 f t 近似为零 下面的列表大致描述了燃烧矿物燃料 特别是汽油 所释放二氧化 碳的源项 f t 其数值在 1850 年后开始变得显著 在本世纪末达 到峰值 然后降低直到燃料消耗完 年使用率年使用率 10000 0210010 5 18500 0212010 0 19501 021508 0 19804 022253 5 20005 023002 0 20508 025000 0 208010 050000 0 上图显示了源项及其对大气和海洋的影响 下半部的三个图线显示 了大气 浅海及深海中碳的变化 两个碱度值由于其基本保持不变 所以没有绘出 开始时 三个区域中的碳维持在平衡状态 所以其 数量在 1850 年之前几乎没有变化 在 1850 t 2500 之间 图的上 半部表示了燃烧矿物燃料所产生的碳进入系统 下半部分表示了系 统的反应 大气是最先受到影响的 在 500 年后将增加四倍多 大 约有一倍的碳进入浅海并最终转入深海区域内 a 重新绘制上图 使用 pchiptx 插入燃料表格 并使用 ode23tx 在默认的容差下求解此微分方程 b 比较三个区域中 5000 年和 1000 年时碳的含量 c 大气中的二氧化碳在何时达到最大值 d 该方程是适度刚性刚性的 这是由不同时间尺度内多种化学反应的 发生造成的 放大图的某些部分 可以发现由于 ode23tx 采用 了小步长所引起的特征锯齿形态 找到这样的一个区域 e 用 MATLAB 的常微分方程求解程序 包括 ode23 ode45 ode113 ode23s 以及 ode15s 尝试采用不同的容差 使用如下代码并报告其计算代价 odeset RelTol 1 e 6 AbsTol 1 e 6 stats on 哪个方法更适合解决此问题 1 Lighthouse Foundation http www lighthouse foundation org 2 J C G Walker Numerical Adventures with Geochemical Cycles Oxford University Press New York 1991 注 1 pchiptx 是 Matlab 函数 pchip 的一个可查看源代码的版本 包 含在 NCM 程序包中 其作用是是实现三次埃米特插值 用于 从函数的一组采样点生成该函数的一个近似的连续表达式 其 用法可参看 插值 pdf 2 刚性问题的定义可参考任何一本计算方法或微分方程数值解的 书 也可参看 常微分方程 pdf 3 类似 ode45 这样的以 ode 开头的函数是 Ma