数值微分
数值微分(numerical differentiation)
根据函数在一些离散点的函数值，推算它在某点的导数或高阶导数的近似值的方法。通常用差商代替微商，或者用一个能够近似代替该函数的较简单的可微函数（如多项式或样条函数等）的相应导数作为能求导数的近似值。例如一些常用的数值微分公式(如两点公式、三点公式等)就是在等距步长情形下用插值多项式的导数作为近似值的。此外，还可以采用待定系数法建立各阶导数的数值微分公式，并且用外推技术来提高所求近似值的精确度。当函数可微性不太好时，利用样条插值进行数值微分要比多项式插值更适宜。如果离散点上的数据有不容忽视的随机误差，应该用曲线拟合代替函数插值，然后用拟合曲线的导数作为所求导数的近似值，这种做法可以起到减少随机误差的作用。数值微分公式还是微分方程数值解法的重要依据。
7.1 数值微分
7.1.1 差商与数值微分
当函数是以离散点列给出时，当函数的表达式过于复杂时，常用数值微分近似计算
的导数。在微积分中，导数表示函数在某点上的瞬时变化率，它是平均变化率的极限；在几何上可解释为曲线的斜率；在物理上可解释为物体变化的速率。
以下是导数的三种定义形式：
（7.1）
在微积分中，用差商的极限定义导数；在数值计算中返璞归真，导数取用差商（平均变化率）作为其近似值。
最简单的计算数值微分的方法是用函数的差商近似函数的导数，即取极限的近似值。下面是与式（7.1）相应的三种差商形式的数值微分公式以及相应的截断误差。
向前差商
用向前差商（平均变化率）近似导数有：
（7.2）
其中的位置在的前面，因此称为向前差商。同理可得向后差商、中心差商的定义。
由泰勒展开
得向前差商的截断误差：
向后差商
用向后差商近似导数有：（7.3）
与计算向前差商的方法类似，由泰勒展开得向后差商的截断误差：
中心差商
用中心差商（平均变化率）近似导数有：
（7.4）
由泰勒展开
得中心差商的截断误差：
差商的几何意义
微积分中的极限定义，表示在处切线的斜率，
即图7.1中直线的斜率；差商表示过和两
点直线的斜率，是一条过的割线。可见数值微分是用近似值内接弦的斜率代替准确值切线的斜率。
图7.1 微商与差商示意图
例7.1给出下列数据，计算，
0.10
－0.02 0.04 0.06 0.8
5.06 5.07 5.065 5.05 5.055
解：（5.07－5.06）/（0.04－0.02）= 0.5
（5.05－5.07）/（0.08－0.04）= －0.5
（5.05－5.055）/（0.08－0.10）= 0.25
（(0.10) －(0.06)）/（0.10－0.06）= 18.75
设定最佳步长
在计算数值导数时，它的误差由截断误差和舍入差两部分组成。用差商或插值公式近似导数产生截断误差，由原始值的数值近似产生舍入误差。在差商计算中，从截断误差的逼
近值的角度看，越小，则误差也越小；但是太小的会带来较大的舍入误差。怎样选择最佳步长，使截断误差与舍入误差之和最小呢？
一般对计算导数的近似公式进行分析可得到误差的表示式，以中心差商为例，截断误差不超过
而舍入误差可用量估计（证明略），其中是函数的原始值的绝对误差限，总误差
为
当时，总误差达到最小值，即（*）
可以看到用误差的表达式确定步长，难度较大，难以实际操作。
通常用事后估计方法选取步长，例如，记为步长等于的差商计算公式，给定误差界，当时，就是合适的步长。
例7.2对函数，取不同的步长计算，观察误差变化规律，从而确定最佳步长。
解：
误差误差
0.10 0.09 3.1630
3.1622
－0.0048
－0.0040
0.05
0.04
3.1590
3.1588
－0.0008
－0.0006
0.08 0.07 0.06 3.1613
3.1607
3.1600
－0.0031
－0.0025
－0.0018
0.03
0.02
0.01
3.1583
3.1575
3.1550
－0.0001
－0.0007
－
0.0032
表中数据显示，当步长从0.10减少到0.03时，数值微分误差的绝对值从0.0048减少到0.0001，而随着的进一步减少，误差的绝对值又有所反弹，表明当步长小于0.03时，舍入误差起了主要作用。
在实际计算中是无法得到误差的准确数值的，这时以最小为标准确定步长，本例中取= 0.04。
7.1.2 插值型数值微分
对于给定的的函数表，建立插值函数，用插值函数的导数近似函数
的导数。
设为上的节点，给定，以为插值点构造插值多项式，以的各阶导数近似的相应阶的导数，即
当时，
（7.5）误差项为：
例7.3给定，并有，计算。解：作过的插值多项式：
将代入得三点端点公式和三点中点公式：
利用泰勒（Taylor）展开进行比较和分析，可得三点公式的截断误差是。
类似地，可得到五点中点公式和五点端点公式：
7.1.3 样条插值数值微分
把离散点按大小排列成，用关系式构造插值点
的样条函数：
当则当时，可用计算导数。
曲线拟合-正文
用连续曲线近似地刻画或比拟平面上离散点组所表示的坐标之间的函数关系。更广泛地说，空间或高维空间中的相应问题亦属此范畴。在数值分析中，曲线拟合就是用解析表达式逼近离散数据，即离散数据的公式化。实践中，离散点组或数据往往是各种物理问题和统计问题有关量的多次观测值或实验值，它们是零散的，不仅不便于处理，而且通常不能确切和充分地体现出其固有的规律。这种缺陷正可由适当的解析表达式来弥补。
数学表述设给定离散数据
(1)
式中x k为自变量x(标量或向量,即一元或多元变量)的取值；y k为因变量y（标量）的相应值。曲线拟合要解决的问题是寻求与(1)的背景规律相适应解析表达式
(2)
使它在某种意义下最佳地逼近或拟合(1),?(x,b)称为拟合模型；为待定参数,当b)仅在?中线性地出现时，称模型为线性的，否则为非线性的。量
称为在x k处拟合的残差或剩余，衡量拟合优度的标准通常有
式中ωk>0为权系数或权重（如无特别指定，一般取为平均权重，即
(k=1,2,…,m),此时无需提到权）。当参数b)使T(b))或Q(b))达到最小时，相应的(2)分别称为在加权切比雪夫意义或加权最小二乘意义下对(1)的拟合，后者在计算上较简便且最为常用。
模型中参数的确定一般的线性模型是以参数b)为系数的广义多项式，即
,(3)
式中g0,g1,…，g n称为基函数。对诸g j的不同选取可构成多种典型的和常用的线性
模型。从函数逼近的观点来看，式(3)还能近似地体现许多非线性模型的性质。
在最小二乘意义下用线性模型(3)拟合离散点组(1)，参数b可通过解方程组(i=0,…,n)来确定，即解关于b0,b1,…,b n的线性代数方程组
(4)
式中(i,j＝0,1,…，n),
方程组(4)通常称为法方程或正规方程，当m>n时一般有惟一解。
至于非线性模型以及非最小二乘原则的情形,参数b)可通过解非线性方程组或最优化计算中的有关方法来确定（见非线性方程组数值解法、最优化）。
模型的选择对于给定的离散数据(1),需恰当地选取一般模型(2)中函数?(x,b))的类别和具体形式,这是拟合效果的基础。若已知(1)的实际背景规律，即因变量y对自变量x的依赖关系已有表达式形式确定的经验公式，则直接取相应的经验公式为拟合模型。反之，可通过对模型(3)中基函数g0,g1,…,g n(个数和种类)的不同选取，分别进行相应的拟合并择其效果佳者。函数g0,g1,…,g n对模型的适应性起着测试的作用，故又称为测试函数。另一种途径是:在模型(3)中纳入个数和种类足够多的测试函数，借助于数理统计方法中的相关性分析和显著性检验，对所包含的测试函数逐个或依次进行筛选以建立较适合的模型（见回归分析）。当然，上述方法还可对拟合的残差（视为新的离散数据）再次进行，以弥补初次拟合的不足。总之，当数据中变量之间的内在联系不明确时，为选择到相适应的模型，一般需要反复地进行拟合试验和分析鉴别。
科学计算—理论、方法 及其基于MATLAB 的程序实现与分析 微分方程（组）数值解法 §1 常微分方程初值问题的数值解法 微分方程（组）是科学研究和工程应用中最常用的数学模型之一。如揭示质点运动规律的Newton 第二定律： ()()()?????'='==0 00022x t x x t x t F dt x d m （１） 和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题： () ()0 0y y t f ay dt dy =+= （２） 的解为： ()()()τττd f e y e t y t t a at ?-+=00 （３） 但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实：
１、绝大多数的常微分方程和常微分方程组得不到（有限形式的）解析解； ２、实际应用中往往只需要知道常微分方程（组）的解在（人们所关心的）某些点处的函数值（可以是满足一定精度要求的近似值）； 如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。 一般的一阶常微分方程（组）的初值问题是指如下的一阶常微分方程（组）的定解问题： ()()0 00,y t y t t t y t F dt dy f =≤≤= （７） 其中 ()()()()???? ?? ? ??=t y t y t y t y n 21 （８） ()()()()???? ?? ? ??=y t f y t f y t f y t F n ,,,,21 （９） 常微分方程（组）的初值问题通常是对一动态过程（动态系统、动力系统）演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。 §１.1 常微分方程（组）的Cauch 问题数值解法概论
偏导数与全导数偏微分与全微分的关系 Document serial number【UU89WT-UU98YT-UU8CB-UUUT-UUT108】
1。偏导数代数意义偏导数是对一个变量求导,另一个变量当做数对x求偏导的话y就看作一个数，描述的是x方向上的变化率对y求偏导的话x就看作一个数，描述的是y方向上的变化率 几何意义对x求偏导是曲面z=f(x,y)在x方向上的切线对y求偏导是曲面z=f(x,y)在x方向上的切线 这里在补充点。就是因为偏导数只能描述x方向或y方向上的变化情况，但是我们要了解各个方向上的情况，所以后面有方向导数的概念。 2。微分偏增量：x增加时f(x,y)增量或y增加时f(x,y) 偏微分：在d e t a x趋进于0时偏增量的线性主要部分d e t a z=f x（x,y）d e t a x+o(d e t a x) 右边等式第一项就是线性主要部分，就叫做在（x，y）点对x的偏微分这个等式也给出了求偏微分的方法，就是用求x的偏导数求偏微分
全增量：x,y都增加时f(x,y)的增量全微分：根号(detax方+detay方）趋于0时，全增量的线性主要部分同样也有求全微分公式，也建立了全微分和偏导数的关系d z=A d x+B d y其中A就是对x求偏导，B就是对y求偏导 希望楼主注意的是导数和微分是两个概念，他们之间的关系就是上面所说的公式。概念上先有导数，再有微分，然后有了导数和微分的关系公式，公式同时也指明了求微分的方法。 3.全导数全导数是在复合函数中的概念，和上面的概念不是一个系统，要分开。u=a(t),v=b(t) z=f[a(t),b(t)] dz/dt 就是全导数，这是复合函数求导中的一种情况，只有这时才有全导数的概念。 d z/d t=(偏z/偏u)(d u/d t)+(偏z/偏v)(d v/d t) 建议楼主在复合函数求导这里好好看看书，这里分为3种情况。1.中间变量一元就是上面的情况，才有全导数的概念。2.中间变量有多元，只能求偏导 3.中间变两有一元也有多元，还是求偏导。 对于你的题能求对x的偏导数，对y的偏导数，z的全微分，不能求全导数
《微分方程数值解法》 【摘要】自然界与工程技术中的很多现象，可以归结为微分方程定解问题。其中，常微分方程求解是微分方程的重要基础内容。但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、Runge —Kutta 方法、Adams 预估校正法以及勒让德谱方法等，通过具体的算例，结合MA TLAB 求解画图，初步给出了一般常微分方程数值解法的求解过程。同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。 【关键词】 常微分方程 数值解法 MA TLAB 误差分析 引言 在我国高校，《微分方程数值解法》作为对数学基础知识要求较高且应用非常广泛的一门课程，不仅 在数学专业，其他的理工科专业的本科及研究生教育中开设这门课程．近四十年来，《微分方程数值解法》不论在理论上还是在方法上都获得了很大的发展．同时，由于微分方程是描述物理、化学和生物现象的数学模型基础，且它的一些最新应用已经扩展到经济、金融预测、图像处理及其他领域 在实际应用中，通过相应的微分方程模型解决具体问题，采用数值方法求得方程的近似解，使具体问题迎刃而解。 2 欧拉法和改进的欧拉法 2.1 欧拉法 2.1.1 欧拉法介绍 首先，我们考虑如下的一阶常微分方程初值问题 ???==0 0)() ,('y x y y x f y （2--1） 事实上，对于更复杂的常微分方程组或者高阶常微分方程，只需要将x 看做向量，（2--1）就成了一个一阶常微分方程组，而高阶常微分方程也可以通过降阶化成一个一阶常微分方程组。 欧拉方法是解常微分方程初值问题最简单最古老的一种数值方法，其基本思路就是把（2--1）中的导数项'y 用差商逼近，从而将一个微分方程转化为一个代数方程，以便求解。 设在[]b a ,中取等距节点h ，因为在节点n x 点上，由（2--1）可得：
第十章 偏微分方程数值解法 偏微分方程问题，其求解十分困难。除少数特殊情况外，绝 大多数情况均难以求出精确解。因此，近似解法就显得更为重要。本章仅介绍求解各类典型偏微分方程定解问题的差分方法。 §1 差分方法的基本概念 1.1 几类偏微分方程的定解问题 椭圆型方程：其最典型、最简单的形式是泊松（Poisson ）方程 ),(22 2 2y x f y u x u u =??+??=? 特别地，当0),(≡y x f 时，即为拉普拉斯（Laplace ）方程，又 称 为调和方程 22 22=??+??=?y u x u u Poisson 方程的第一边值问题为 ?? ???Ω ?=Γ=Ω∈=??+??Γ∈),(),(),() ,(),(22 22y x y x u y x y x f y u x u y x ?
其中 Ω为以Γ为边界的有界区域，Γ为分段光滑曲线，ΓΩY 称为定解区域，),(y x f ，),(y x ?分别为Ω，Γ上的已知连 续函数。 第二类和第三类边界条件可统一表示为 ),(),(y x u u y x ?α=??? ? ??+??Γ∈n 其中n 为边界Γ的外法线方向。当0=α时为第二类边界条件， 0≠α时为第三类边界条件。 抛物型方程：其最简单的形式为一维热传导方程 2 20(0)u u a a t x ??-=>?? 方程可以有两种不同类型的定解问题： 初值问题 ?? ???+∞ <<∞-=+∞<<-∞>=??-??x x x u x t x u a t u )()0,(,00 22? 初边值问题
2 212 00,0(,0)()0(0,)(),(,)()0u u a t T x l t x u x x x l u t g t u l t g t t T ????-=<<<
郑州大学研究生课程（2012-2013学年第一学期）数值分析 Numerical Analysis 习题课 第八章常微分方程数值解法
待求解的问题：一阶常微分方程的初值问题/* Initial-Value Problem */: ?????=∈=0 )(] ,[),(y a y b a x y x f dx dy 解的存在唯一性（“常微分方程”理论）：只要f (x , y ) 在[a , b ] ×R 1 上连续，且关于y 满足Lipschitz 条件，即存在与x , y 无关的常数L 使 对任意定义在[a , b ] 上的y 1(x ) 和y 2(x ) 都成立，则上述IVP 存在唯一解。 1212|(,)(,)||| f x y f x y L y y ?≤?一、要点回顾
§8.2 欧拉(Euler)法 通常取(常数),则Euler 法的计算格式 h h x x i i i ==?+1?? ?=+=+) (),(001x y y y x hf y y i i i i i =0,1,…,n ( 8.2 )
§8.2 欧拉(Euler)法(1) 用差商近似导数 )) (,()()()()(1n n n n n n x y x hf x y x y h x y x y +=′+≈+?? ?=+=+) (),(01a y y y x hf y y n n n n 差分方程初值问题向前Euler 方法h x y x y x y n n n ) ()()(1?≈ ′+)) (,() ()(1n n n n x y x f h x y x y ≈?+))(,()(n n n x y x f x y =′
1。偏导数 代数意义 偏导数是对一个变量求导,另一个变量当做数 对x求偏导的话y就看作一个数，描述的是x方向上的变化率 对y求偏导的话x就看作一个数，描述的是y方向上的变化率 几何意义 对x求偏导是曲面z=f(x,y)在x方向上的切线 对y求偏导是曲面z=f(x,y)在x方向上的切线 这里在补充点。就是因为偏导数只能描述x方向或y方向上的变化情况，但是我们要了解各个方向上的情况，所以后面有方向导数的概念。 2。微分 偏增量：x增加时f(x,y)增量或y增加时f(x,y) 偏微分：在detax趋进于0时偏增量的线性主要部分 detaz=fx（x,y）detax+o(detax) 右边等式第一项就是线性主要部分，就叫做在（x，y）点对x的偏微分 这个等式也给出了求偏微分的方法，就是用求x的偏导数求偏微分 全增量：x,y都增加时f(x,y)的增量 全微分：根号(detax方+detay方）趋于0时，全增量的线性主要部分
同样也有求全微分公式，也建立了全微分和偏导数的关系 dz=Adx+Bdy 其中A就是对x求偏导，B就是对y求偏导 希望楼主注意的是导数和微分是两个概念，他们之间的关系就是上面所说的公式。概念上先有导数，再有微分，然后有了导数和微分的关系公式，公式同时也指明了求微分的方法。 3.全导数 全导数是在复合函数中的概念，和上面的概念不是一个系统，要分开。 u=a(t),v=b(t) z=f[a(t),b(t)] dz/dt 就是全导数，这是复合函数求导中的一种情况，只有这时才有全导数的概念。 dz/dt=(偏z/偏u)(du/dt)+(偏z/偏v)(dv/dt) 建议楼主在复合函数求导这里好好看看书，这里分为3种情况。1.中间变量一元就是上面的情况，才有全导数的概念。2.中间变量有多元，只能求偏导 3.中间变两有一元也有多元，还是求偏导。 对于你的题能求对x的偏导数，对y的偏导数，z的全微分，不能求全导数 如果z=f(x^2,2^x) 只有这种情况下dz/dx才是全导数！
微分方程的分类及其数值解法 微分方程的分类： 含有未知函数的导数，如dy/dx=2x 、ds/dt=0.4都是微分方程。 一般的凡是表示未知函数、未知函数的导数与自变量之间的关系的方程，叫做微分方程。未知函数是一元函数的，叫常微分方程；未知函数是多元函数的叫做偏微分方程。微分方程有时也简称方程。 一、常微分方程的数值解法： 1、Euler 法： 00d (,), (1.1)d (), (1.2) y f x y x y x y ?=???=? 001 (),(,),0,1,,1n n n n y y x y y hf x y n N +=??=+=-? (1.4) 其中0,n b a x x nh h N -=+=. 用(1.4)求解(1.1)的方法称为Euler 方法。 后退Euler 公式???+==+++),,(),(111 00n n n n y x hf y y x y y 梯形方法公式 )].,(),([2 111+++++=n n n n n n y x f y x f h y y 改进的Euler 方法11(,),(,),1().2p n n n c n n p n p c y y hf x y y y hf x y y y y ++?=+??=+???=+??? 2、Runge-Kutta 方法： p 阶方法 ： 1()O h -=?总体截断误差局部截断误差 二阶Runge-Kutta 方法 ??? ????++==++=+),,(),,(,2212 1211hk y h x f k y x f k k h k h y y n n n n n n
图５畅２ 令珔h ＝h λ，则y n ＋１＝１＋珔 h ＋１２珔h ２ ＋１６珔h ３＋１２４ 珔 h ４y n ．由此可知，绝对稳定性区域在珔h ＝h λ复平面上满足 ｜１＋珔 h ＋１２珔h ２＋１６珔h ３＋１２４珔h ４ ｜≤１的区域，也就是由曲线 １＋珔h ＋ １２珔h ２＋１６珔h ３＋１２４ 珔h ４＝e i θ 所围成的区域．如图５畅２所示． 例22 用Euler 法求解 y ′＝－５y ＋x ，y （x ０）＝y ０， x ０≤x ≤X ． 从绝对稳定性考虑，对步长h 有何限制？ 解 对于模型方程y ′＝λy （λ＜０为实数）这里λ＝抄f 抄y ＝－５．由 ｜１＋h λ｜＝｜１－５h ｜＜１ 得到对h 的限制为：０＜h ＜０畅４． 四、习题 １畅取步长h ＝０畅２，用Euler 法解初值问题 y ′＝－y －x y ２ ， y （０）＝ １． （０≤x ≤０畅６）， ２畅用梯形公式解初值问题 y ′＝８－３y ， （１≤x ≤２），
取步长h＝０畅２，小数点后至少保留５位． ３畅用改进的Euler公式计算初值问题 y′＝１x y－１x y２， y（１）＝０畅５， １＜x＜１畅５， 取步长h＝０畅１，并与精确解y（x）＝ x １＋x比较． ４畅写出用梯形格式的迭代算法求解初值问题 y′＋y＝０， y（０）＝１ 的计算公式，取步长h＝０畅１，并求y（０畅２）的近似值，要求迭代误差不超过１０－５． ５畅写出用四阶经典Runge唱Kutta法求解初值问题 y′＝８－３y， y（０）＝２ 的计算公式，取步长h＝０畅２，并计算y（０畅４）的近似值，小数点后至少保留４位． ６畅证明公式 y n＋１＝y n＋h９（２K１＋３K２＋４K３）． K１＝f（x n，y n）， K２＝f x n＋h２，y n＋h２K１， K３＝f x n＋３４h，y n＋３４h K２， 至少是三阶方法． ７畅试构造形如 y n＋１＝α（y n＋y n－１）＋h（β０f n＋β１f n－１）
并与真解u(x) 2e x x 1相比较. 微分方程数值方法 常微分方程初值问题习题一 u' ax b, u(0) 0, 分别写出Euler 法和改进的Euler 法的近似解 府 的表达式，并求 它们与真解u(x) -ax 2 bx 的差u(X m ) U m . 2. 取步长h 0.1，分别用Euler 法和改进的Euler 法求下列初值问 题的解，并与真解相比较. 真解 u(x) .1 2x ； 2 ，u x . c (2) u 2 ,1 x 2, x u u(1) 2, 1 真解 u(x) x(8 31 n x)3 ； u x u '広乔 u(1) 1, 3 1 真解 u(x) (4x 2 3x 2)3. X 2 3. 用Euler 法计算0￡dt 在x 0.1，0.2的近似值. 4. 取步长h 0.2，用四阶Runge-Kutta 法解 u' u x, 0 x 1, u(0) 1, 1.对初值问题 (1) u' u 2x 0x1, u(0) u 1 , (3) 1 x 1.5,
5. 设 f(x,u)关于 u 满足 Lipschitz 条件，证明 N 级 Runge-Kutta 法中的增量函数 (x,u,h)关于u 也满足 Lipschitz 条件. 6. 对初值问题 u' u x 1, u(0) 1, 写出四阶Taylor 级数法和四阶 Runge-Kutta 法的计算公式，它们 是否相同. 7. 证明改进的Euler 法的绝对稳定区间是(-2，0). 8.证明：当h( h)满足 3 4 h h 24 时，四阶 Runge-Kutta 法绝对稳定. 9. 用Tayor 展开确定下面多步法中的系数，使其阶尽可能高，并求 局部截断误差的主项. 10. 对初值问题 u'' f(x,u), u(X °) u °,u'(x 0) u 10, 确定求解公式 (3) u m1 a 2u m 1 h( m 1 2 ). (1) u m 1 a 1u m a 2u m 1 h 0 f m 1 ；
常微分方程数值方法 1、欧拉方法：1,,1,0),,(1-=+=+n k y t hf y y k k k k . function E=euler(f,a,b,ya,n) % Input - f is the function entered as a string 'f' % - a and b are the left and right end points % - ya is the initial condition y(a) % - n is the number of steps % Output - E=[T' Y'] where T is the vector of abscissas and % Y is the vector of ordinates h=(b-a)/n; T=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; Y(1)=ya; for j=1:n Y(j+1)=Y(j)+h*feval(f,T(j),Y(j)); end E=[T' Y']; 【例】 用欧拉方法求解区间]3,0[内的初值问题：1)0(,2'=-=y y t y 。 f=inline('(t-y)/2','t','y');a=0;b=3;ya=1;n=12; %n=3,6,12,24,48,96... E=euler(f,a,b,ya,n),plot(E(:,1),E(:,2),'r*'),hold on 符号解：y=dsolve('Dy=(t-y)/2','y(0)=1') h=(3-0)/12;t=0:h:3;y=eval(y);[t' y'] 用图比较数值解：(f 为ode 函数文件) ode45('f',[0,3],1) 2、休恩(Huen)方法(即改进Euler 方法)： 1 ,,1,0)],,(,(),([211-=+++=++n k y t hf y t f y t f h y y k k k k k k k k function H=heun(f,a,b,ya,n) % Input - f is the function entered as a string 'f' % - a and b are the left and right end points % - ya is the initial condition y(a) % - n is the number of steps % Output - H=[T' Y'] where T is the vector of abscissas and % Y is the vector of ordinates h=(b-a)/n; T=zeros(1,n+1); Y=zeros(1,n+1); T=a:h:b; Y(1)=ya; for j=1:n k1=feval(f,T(j),Y(j)); k2=feval(f,T(j+1),Y(j)+h*k1); Y(j+1)=Y(j)+(h/2)*(k1+k2); end H=[T' Y'];
偏导数与全微分习题 1. 设y x y x y x f arcsin )1(),(-+=，求)1,(x f x '。 2. 习题8 17题。 3. 设?? ??? =+≠++=0 001sin ),(22222 2 y x y x y x y y x f ，考察f (x , y )在点（0,0）的偏导数。 4. 考察?? ??? =+≠++=0 001sin ),(22222 2 y x y x y x xy y x f 在点 （0,0）处的可微性。 5. 证 明 函 数 ?? ???=+≠+++=0 001sin )(),(222 22 22 2y x y x y x y x y x f 在 点（0,0）连续且偏导数存在，但偏导数在（0,0）不连续，而f (x , y )在点（0,0）可微。 }
1． 设y x y x y x f arcsin )1(),(-+=，求)1,(x f x '。 y y x y x y y x f x 1) (2111 )1(1),(21 ??- -+='- ∴ 1)1,(='x f x 。 ： &
2.习题8 17题。 17. 设22)()(ln b y a x z -+-=（a , b 为常数），证明 02 22 2=??+??y z x z 。 先化简函数 ))()ln((2 1 22b y a x z -+-=， ， 2 222)()() ()()()(221b y a x a x b y a x a x x z -+--= -+--?=??， 2222) ()() ()()()(221b y a x b y b y a x b y y z -+--=-+--?=??， 2 22 2 222 2))()(()(2)()(b y a x a x b y a x x z -+----+-= ?? 2 22 22) )()(()()(b y a x a x b y -+----= ， 2 222 222 2))()(()(2)()(b y a x b y b y a x y z -+----+-= ?? 2 2222) )()(()()(b y a x b y a x -+----= ， ∴ 02 22 2=??+ ??y z x z 。 3. $
i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法，其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程，再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<< <<= （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。
《偏微分方程数值解法》 课程设计 题目: 六点对称差分格式解热传导方程的初边 值问题 姓名: 王晓霜 学院: 理学院 专业: 信息与计算科学 班级: ０91１01２ 学号: ０９１10１218 指导老师：翟方曼 ２０12年12月14
日 一、题目 用六点对称差分格式计算如下热传导方程的初边值问题 222122,01,01(,0),01 (0,),(1,),01x t t u u x t t x u x e x u t e u t e t +???=<<<≤?????=≤≤??==≤≤??? 已知其精确解为 2(,)x t u x t e += 二、理论 1.考虑的问题 考虑一维模型热传导方程 （1.１) )(22x f x u a t u +??=??,T t ≤<0 其中a 为常数。)(x f 是给定的连续函数。（1.１)的定解问题分两类: 第一,初值问题（Cauch ｙ 问题)：求足够光滑的函数()t x u ,,满足方程（１.1）和初始条件： (1.２) ()()x x u ?=0,， ∞<<∞-x 第二,初边值问题（也称混合问题）:求足够光滑的函数()t x u ,，满足方程(１.１)和初始条件： ()13.1 ()()x x u ?=0,， l x l <<- 及边值条件 ()23.1 ()()0,,0==t l u t u ， T t ≤≤0 假定()x f 和()x ?在相应的区域光滑,并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。 现在考虑边值问题（１.1),（1.３)的差分逼近 取 N l h = 为空间步长,M T =τ为时间步长,其中N ,M 是自然数,
x 1 ?若步长趋于零时，差分方程的截断误差 R m 0,则差分方程的解 U i m 趋近于微分方 程的解U m ?此结论 ________ (错或对)； 1 2.一 阶 Sobolev 空间 H ( ) f (x,y) f , f x , f y L ?() 关于内积(f,g )1 _____________________________________ 是Hilbert 空间； 3 ?对非线性(变系数)差分格式，常用 ____________ 系数法讨论差分格式的 ________ 稳定性； 4?写出y x 3在区间［1,2］上的两个一阶广义导数： ______________________________________ _____ ____ ______________ _ ____ ________ ; 5 ?隐式差分格式关于初值是无条件稳定的 ?此结论 ________ (错或对)。 (13分)设有椭圆型方程边值问题 0.1作正方形网格剖分 。 (1) 用五点菱形差分格式将微分方程在内点离散化； (2) 用截断误差为 O (h 2)的差分法将第三边界条件离散化; (3) 整理后的差分方程组为 U C 三.(12)给定初值问题 u x,0 x 1 取时间步长 0.1，空间步长h 0.2。试合理选用一阶偏心差分格式(最简显格式) 2 u ~2 x 2 u ~2 y 0 x 0.3 0.2 x 0.3 2y 1, — u n 2x y 0.2
并以此格式求出解函数u(x,t)在x 0.2,t 0.2处的近似值。 x
1.所选用的差分格式是： 2 .计算所求近似值： 1 a k 1 四.（12分）试讨论差分方程 u l 1 k k k 1 u | r u | 1 u | , r h a 1 h 逼近微分方程 u a u 0 t x 的截断误差阶R 。 思路一：将r 带入到原式，展开后可得格式是在点（ l+1/2,k+1/2 ）展开的。 思路二：差分格式的用到的四个点刚好是矩形区域的四个顶点，可由此构造中心点的差分格 式。 2 —2 ,考虑 Du Fort-Frankel 格式 X 试论证该格式是否总满足稳定性的 Von-Neumann 条件？ 六. （12分）（1 ）由Green 第一公式推导 Green 第二公式: （2） 对双调和方程边值问题 n 2 选择函数集合（空间）为： 推导相应的双线性泛函和线性泛函: A （u,v ） F （v ） 相应的虚功问题为： 极小位能问题为 七. （ 12分）设有常微分方程边值问题 y y f （x ） , a x b y a 1, y b 1 五.（12分） 对抛物型方程 U |k1 U |k 2 |k 1 (U |k1 U |k1) U |k 1 ) 2 (u)vdxdy G (u) u vdxdy :[v v u ]ds n f (x,y) (x,y) g 1(x , y), g 2(x, y) (x,y),
微分方程数值解法
《微分方程数值解法》论文 常微分方程初值问题数值解法的讨论与比较一、一阶常微分方程的初值问题 科学技术中常常要求解常微分方程的定解问题，这类问题最简单的形式是一阶方程的初值问题
?Skip Record If...? 我们知道，只要函数适当光滑，譬如关于满足Lipschitz条件 （1.3） 理论上就可以保证初值问题(1.1),(1.2)的解存在并且唯一。 虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法。 二、欧拉法 我们知道，在平面上，微分方程(1.1)的解称作它的积分曲线。积分曲线上一点的切线斜率等于函数的值，如果按函数在平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致，基于上述几何解释，从初始点出发，先依方向场在该点的方向推进到上一点，然后再从依方向场的方向推进到上一点，循此前进推出一条折线，一般地，设已做出该折线的顶点，过依方向场的方向再推进到，显然两个顶点，的坐标有关系 即 (2.1)
这就是著名的欧拉(Euler)公式。若初值已知，则依公式(2.1)可逐步算出 例1求解初值问题 ?Skip Record If...? 解欧拉公式的具体形式为 取步长，计算结果如下表： 表1 计算结果对比 0.1 0.2 0.3 0.4 0.5 1.1000 1.1918 1.2774 1.3582 1.4351 1.0954 1.1832 1.2649 1.3416 1.4142 0.6 0.7 0.8 0.9 1.0 1.5090 1.5803 1.6498 1.7178 1.7848 1.4832 1.5492 1.6125 1.6733 1.7321 初值问题( 2.2)有解，按这个解析式子算出的准确值同近似值一起列在表1，两者相比较可以看出欧拉方法的精度很差。
1. 课本2p 有证明 2. 课本812,p p 有说明 3. 课本1520,p p 有说明 4. Rit2法，设n u 是u 的n 维子空间，12,...n ???是n u 的一组基底，n u 中的任一元素n u 可 表为1n n i i i u c ?==∑ ，则,11 11()(,)(,)(,)(,)22j n n n n n n i j i j j i j j J u a u u f u a c c c f ???=== -=-∑∑是12,...n c c c 的二次函数，(,)(,)i j j i a a ????=，令 () 0n j J u c ?=?，从而得到12,...n c c c 满足1 (,)(,),1,2...n i j i j i a c f j n ?? ?===∑，通过解线性方程组，求的i c ，代入1 n n i i i u c ?==∑， 从而得到近似解n u 的过程称为Rit2法 简而言之，Rit2法：为得到偏微分方程的有穷维解，构造了一个近似解，1 n n i i i u c ? == ∑， 利用,11 11()(,)(,)(,)(,)22j n n n n n n i j i j j i j j J u a u u f u a c c c f ???===-=-∑∑确定i c ，求得近似解n u 的过程 Galerkin 法：为求得1 n n i i i u c ? == ∑形式的近似解，在系数i c 使n u 关于n V u ∈，满足(,)(,) n a u V f V =，对任 意 n V u ∈或（取 ,1j V j n ?=≤≤） 1 (,)(,),1,2...n i j i j i a c f j n ???===∑的情况下确定i c ，从而得到近似解1 n n i i i u c ?==∑的过程称 Galerkin 法为 Rit2-Galerkin 法方程： 1 (,)(,)n i j i j i a c f ???==∑ 5. 有限元法：将偏微分方程转化为变分形式，选定单元的形状，对求解域作剖分，进而构 造基函数或单元形状函数，形成有限元空间，将偏微分方程转化成了有限元方程，利用 有效的有限元方程的解法，给出偏微分方程近似解的过程称为有限元法。 6. 解：对求解区间进行网格剖分，节点01......i n a x x x x b =<<<<=得到相邻节点1,i i x x -
第十四章 偏导数与全微分 §1. 偏导数与全微分的概念 1．求下列函数的偏导数： (1) 2 2 2 ln()u x x y =+； (2) ()cos()u x y xy =+； (3) arctan x u y =； (4) sin()xy u xye =. 2．设22 22 221sin , 0,(,)0, 0.y x y x y f x y x y ?+≠?+=??+=? ，考察函数在(0,0)点的偏导数. 3 ．证明函数u =(0,0)点连续但偏导数不存在. 4．求下列函数的全微分： (1) u = (2) yz x u xe e y -=++.
5．求下列函数在给定点的全微分： (1) u =在点(1,1,1)； (2) (u x y =+-0，1）. 6．证明函数22222 22, 0,(,) 0, 0.x y x y f x y x y x y ?+≠?=+??+=? 在(0,0)点连续且偏导数存在，但在此点不可微。 7 ．证明：函数22 220(,)0, 0x y f x y x y +≠=+=?在点(0, 0)处偏导数存在，但不可微. 8．设,x y 很小，利用全微分推出下列式(1)(1)m n x y ++的近似公式：
9．求下列函数指定阶的偏导数： (1) 3 3 sin sin u x y y x =+,求633u x y ???； (2) ln()u ax by =+,求m n m n u x y +???. §2. 求复合函数偏导数的链式法则 1．求下列函数指定的偏导数： （1）．设(,,),x y z Φ=Φ ,,,x u v y u v z uv =+=-=求, u v ?Φ?Φ ??. （2） 设),,22(xyz z y x f z --=求x z ?? 2. 求下列函数指定的偏导数（假定所有二阶偏导数都连续） (1) 2 2 (,)u f xy x y =，22u x ?? ； (2) (,)x y u f y z =,2u x y ???； (3) 2 2 2 ()u f x y z =++,22u y ??； (4) (,,)x u f x y xy y =+,2u y x ???.
包括基本概念，差分格式的构造、截断误差和稳定性，这些内容是贯穿整个教材的主线。解答问题关键在过程，能够显示出你已经掌握了书上的内容，知道了解题方法。这次考试题目的类型：20分的选择题，主要是基本概念的理解，后面有五个大题，包括差分格式的构造、截断误差和稳定性。 习题一 1． 略 2． y y x f -=),(，梯形公式：n n n n n n y h h y y y h y y )121(),(2111+-+=+- =+++，所以0122)1(01])121[()121()121(y h h y h h y h h y h h n h h n n n +--+--+-+=+-+==+-+= ，当0→h 时， x n e y -→。 同理可以证明预报-校正法收敛到微分方程的解. 3． 局部截断误差的推导同欧拉公式; 整体截断误差: ? ++++++-++≤1 ),())(,(11111n n x x n n n n n n n dx y x f x y x f R εε 11)(++-++≤n n n y x y Lh R ε，这里R R n ≤ 而111)(+++-=n n n y x y ε，所以 R Lh n n += -+εε1)1(，不妨设1 Matlab 应用 使用Euler 和Rungkutta 方法解臂状摆的能量方程 背景 单摆是常见的物理模型，为了得到摆角θ的关于时间的函数，来描述单摆运动。由角动量定理我们知道 ε J M = 化简得到 0sin 22 =+θθl g dt d 在一般的应用和计算中，只考虑摆角在5度以内的小摆动，因为可以吧简化为θ，这样比较容易解。实际上这是一个解二阶常微分方程的问题。 在这里的单摆是一种特别的单摆，具有均匀的质量M 分布在长为2的臂状摆上， 使用能量法建立方程 W T = h mg ?=2J 2 1ω 化简得到 θθcos 35499.722=dt d 重力加速度取9.80665 1使用欧拉法 令dx dy z =，这样降阶就把二阶常微分方程转化为一阶微分方程组，再利用向前Euler 方法数值求解。 y(i+1)=y(i)+h*z(i); z(i+1)=z(i)+h*7.35499*cos(y(i)); y(0)=0 z(0)=0 精度随着h 的减小而更高，因为向前欧拉方法的整体截断误差与h 同阶，（因为是用了泰勒公式）所以欧拉方法的稳定区域并不大。 2.RK4-四阶龙格库塔方法 使用四级四阶经典显式Rungkutta公式 稳定性很好，RK4法是四阶方法，每步的误差是h5阶，而总积累误差为h4阶。所以比欧 拉稳定。 运行第三个程序：在一幅图中显示欧拉法和RK4法，随着截断误差的积累，欧拉法产生了 较大的误差 h=0.01 h=0.0001，仍然是开始较为稳定，逐渐误差变大 总结：RK4是很好的方法，很稳定，而且四阶是很常用的方法，因为到五阶的时候精度并没有相应提升。通过这两种方法计算出角度峰值y=3.141593，周期是1.777510。 第十三讲：多元函数的偏导数与全微分的练习题答案 一、单项选择题（每小题4分，共24分） 1． 设2(,)f x y x y xy y +-=+ 则(,)f x y = （A ） A ． ()2x x y - B ．2xy y + C ．()2 x x y + D ．2x xy - 解： (,)()f x y x y x y y +-=+ []1()()()2 x y x y x y = ++-- (,)()2x f x y x y ∴=- 2． 22 1cos lim 1x x y o e y x y →→++= （D ） A ． 0 B ．1 C ． 1e D ． 2 e 解：22cos (,)1x e y f x y x y =++在点（1，0）连续 '221cos cos 0lim 11102x x y o e y e e x y →→∴==++++ 3．设(,) f x y 在点00(,)x y 处有偏导数存在，则0000(2,)(,)lim h o f x h y f x h y h →+--=（D ） A ．0 B ．'00(,)x f x y C ．'002(,)x f x y D ．'003(,)x f x y 解：原式=0000(2,)(,)lim 22h o f x h y f x y h →+-? 0000(,)(,)lim h o f x h y f x y h →--+- ='''0000002(,)(,)3(,)x x x f x y f x y f x y += 4．(,)z f x y =偏导数存在是(,)z f x y =可微的 （B ） A ．充分条件 B ．必要条件 C ．充分必要条件 D ．无关条件数值分析Matlab作业龙格库塔欧拉方法解二阶微分方程
第十三讲：多元函数的偏导数与全微分的练习题答案
相关主题
文本预览