常微分方程
常微分方程
目录第一章理论部分1常微分方程数值解简介1.1欧拉方法1.2龙格-库塔方法?2Matlab简介62.1Matlab的主要特点2.2Matlab操作和桌面简介2.3Matlab图形功能.2.4Matlab的M文件设计.14第二章实验部分·181欧拉方法及其Matlab程序..1.1显式欧拉方法的Matlab程序181.2隐式欧拉方法的Matlab程序-181.3改进欧拉方法的Matlab程序·19232龙格-库塔方法及其Matlab程序.232.1二阶龙格--库塔方法及其Matlab程序2.2四阶龙格--库塔方法及其Matlab程序24
目 录 第一章 理论部分 1 常微分方程数值解简介 . 1 1.1 欧拉方法 . 1 1.2 龙格-库塔方法 . 3 2 Matlab 简介 . 6 2.1 Matlab 的主要特点 .6 2.2 Matlab 操作和桌面简介 .6 2.3 Matlab 图形功能 .9 2.4 Matlab 的 M 文件设计.14 第二章 实验部分 1 欧拉方法及其 Matlab 程序.18 1.1 显式欧拉方法的 Matlab 程序 .18 1.2 隐式欧拉方法的 Matlab 程序 .18 1.3 改进欧拉方法的Matlab程序 .19 2 龙格-库塔方法及其 Matlab 程序 .23 2.1 二阶龙格-库塔方法及其 Matlab 程序 .23 2.2 四阶龙格-库塔方法及其 Matlab 程序 .24
前言《常微分方程》是数学专业的学生必修的一门专业基础核心课，也是自然科学和工程技术各领域中应用广泛的数学工具。在常微分方程的教学中，主要集中于求解常微分方程的解析解，而在实际应用中有重大意义的许多微分方程，即使它们能满足非常广泛的解的存在唯一性条件，但它们的解常常不能表达成初等函数的形式即解析解。而对这类方程实际上也不需要精确的解析解，有时只需一些从初值出发的特殊点上的近似值，对这类问题的讨论最常用的方法就是数值积分，即对微分方程进行数值解。目前，各类教材中编写的数学软件应用主要集中于求微分方程的解析解，而对微分方程的数值解没有系统的阐述，出于此目的我们编写了本教材。本教材按照常微分方程实验课教学大纲的要求，在编写过程中，力求通过实验培养学生使用数学软件解决微分方程数值解的能力，使学生掌握软件的命令、功能、计算，把数学知识和科学计算软件有机的结合起来，给出问题的解决方案。本教材共有四章，分为两大部分：第一部分概括的讲述了微分方程的数值解的基本概念及Matlab数学软件。第二部分系统的讲述了Euler方法求微分方程的数值解以及Runge-Kutta方法求微分方程的数值解的Matlab程序。由于时间比较仓促，加之编者水平有限，在内容的取舍和安排上考虑不周之处一定会有，甚至缺点错误也在所难免，恩请同行和读者不各斧正。作者在撰写过程中，得到了内蒙古农业大学理学院领导和数学系领导的支持和帮助，作者在此致以衷心的谢意！编者2018年7月
前言 《常微分方程》是数学专业的学生必修的一门专业基础核心课，也是自然 科学和工程技术各领域中应用广泛的数学工具。在常微分方程的教学中，主要集 中于求解常微分方程的解析解，而在实际应用中有重大意义的许多微分方程，即 使它们能满足非常广泛的解的存在唯一性条件，但它们的解常常不能表达成初等 函数的形式即解析解。而对这类方程实际上也不需要精确的解析解，有时只需一 些从初值出发的特殊点上的近似值，对这类问题的讨论最常用的方法就是数值积 分，即对微分方程进行数值解。目前，各类教材中编写的数学软件应用主要集中 于求微分方程的解析解，而对微分方程的数值解没有系统的阐述，出于此目的， 我们编写了本教材。 本教材按照常微分方程实验课教学大纲的要求，在编写过程中，力求通过实 验培养学生使用数学软件解决微分方程数值解的能力，使学生掌握软件的命令、 功能、计算，把数学知识和科学计算软件有机的结合起来，给出问题的解决方案。 本教材共有四章，分为两大部分：第一部分概括的讲述了微分方程的数值解 的基本概念及 Matlab 数学软件。第二部分系统的讲述了 Euler 方法求微分方程的 数值解以及 Runge-Kutta 方法求微分方程的数值解的 Matlab 程序。 由于时间比较仓促，加之编者水平有限，在内容的取舍和安排上考虑不周之 处一定会有，甚至缺点错误也在所难免，恳请同行和读者不吝斧正。 作者在撰写过程中，得到了内蒙古农业大学理学院领导和数学系领导的支持 和帮助，作者在此致以衷心的谢意！ 编者 2018 年 7 月
第一章 理论部分1常微分方程数值解简介[=f(x,y)初值问题(1）的数值解法，是通过微分方程离散化而给dx[y(a)=α出解在某些节点上的近似值。在[a,b]上引入节点()a=<x..<x,=b，h=x--(k=],,n)称为步长。在多数情况下，采用等步长，即h=-二,x=α+h(k=0.1,m)。n记初值问题的为准确解为(x)，记(x）的近似值为yk，记f(xy）为f。1.1欧拉方法1.1.1显式欧拉方法设节点为α=x。<x<x，=b，初值问题（1）的显式欧拉方法为[ = a[yk+=yk+hk,=0,],..,n-1其中h=X+-X,=f(xk)。1.1.2隐式欧拉方法和梯形方法若将(）在x+1展开y"(),x≤n≤Xk+y(x)= y(xk+1)-hf(Xk+1,(xk+1)+21忽略h?项，用ykyk+和fk+=f(xk1,yk+）分别近似(x），(xk+）及f(xk+1,(xk+）)，可以得另一计算公式Yk1 = yx +h,f(xk+1,k+1),k=0,l,-.-,n-1称为隐式Enler方法。在显式和稳式Euler方法中，忽略的项都是h?项，为了得到更高精确度的方法，我们可将1y(5)h,X≤5k≤Xkly(Xk+1)= y(x)+ hef(xk+1,J(x))+ 21
1 第一章 理论部分 1 常微分方程数值解简介 初值问题 (, ) (1) ( ) dy f xy dx y a α ⎧ ⎪ = ⎨ ⎪ ⎩ = 的数值解法，是通过微分方程离散化而给 出解在某些节点上的近似值。 在[ ] a,b 上引入节点{ } 0 0 1 : n k n k x ax x x b = = <<< = " ， 1( 1, , ) k kk hxxk n = − = − " 称为步长。在多数情况下，采用等步长，即 , x a kh(k 0,1 ,n) n b a h k = + = " − = 。 记初值问题的为准确解为 y(x),记 ( ) k y x 的近似值为 k y ，记 ( , ) k k f x y 为 k f 。 1.1 欧拉方法 1.1.1 显式欧拉方法 设节点为a = x0 < x1" < xn = b ，初值问题（1）的显式欧拉方法为 ⎩ ⎨ ⎧ = + = − = + , 0,1, , 1 1 0 y y h f k n y a k k k k " 其中 , ( , ) k k 1 k k k k h = x − x f = f x y + 。 1.1.2 隐式欧拉方法和梯形方法 若将 ( ) k y x 在 xk +1展开 1 1 () ( ) ( k k kk y x yx hf x = − + + , 2 1 1 ( )) ( ) , 2! k kk yx y h + + ′′ η k ≤ k ≤ k+1 x η x 忽 略 2 h 项，用 1 , k k+ y y 和 ( , ) k+1 = k+1 k+1 f f x y 分别近似 ( ) k y x ， ( ) k+1 y x 及 ( , ( )) k+1 k+1 f x y x ，可以得另一计算公式 yk +1 = yk + hk f (xk +1 , yk+1 ), k = 0,1,",n −1 称为隐式 Enler 方法。 在显式和稳式 Euler 方法中，忽略的项都是 2 h 项，为了得到更高精确度的方 法，我们可将 1 2 1 1 ( ) , 2 1 ( ) ( ) ( , ( )) + + ≤ ≤ + = + + ′′ k k k k k k k k k k y x y x h f x y x y ξ h x ξ x
(x+1)= y(x)+hf(xk+),y(x+)-y"()h,≤n≤xk+1取平均得+()+()+)(y(Xk+1)= y(x) +当(x)三次连续可微时，y(）-y"(n）=O(he）。忽略O(h）项，用yk,k分+(x)+()）称为梯形方法。别近似y(x),(xk+1)，得yk+1=yk+21.1.3改进欧拉方法在实际计算中，先用显式Euler方法所得的k为y()，再用梯形方法改进一次[Jk+ = y +hf(x,yk)[f(x,yx)+(x++,J+)Yk+I=yk +k=0,1,,n-1称为改进Euler方法，改进Euler方法还可写成：+()+(+()Yk+I = Yk2hkk,+hkk或k, =f(xk,yk)k, =f(xk+,y+hk,)Yk+l=y+02M21.1.4单步法的局部截断误差、整体截断误差设所用数值法为[y=α[yk+=yk+hd(x,yk,yk+1,h,),k=0,1..",n-1设(x)是(1)的准确解,称Tk+=(xk+1)-(x）-h(xk,(x),(xk+1),h）为此数值法在xk+的局部截断误差；称e=(x）-yk(k=0,l,,n)为此数值法在x点的整体截断误差;如果对充分小的h>0,x,=a+kh,成立maxly(x)-y|≤Chp(p≥1)常数C独立于h，就称此方法是p阶方法2
2 1 2 1 1 1 ( ) , 2 1 ( ) ( ) ( , ( )) + + + ≤ ≤ + = + − ′′ k k k k k k k k k k y x y x h f x y x y η h x η x 取平均得 [ ( ) ( )] 4 [ ( , ( )) ( , ( ))] 2 ( ) ( ) 2 1 1 1 k k k k k k k k k k y y h f x y x f x y x h y x + = y x + + + + + ′′ ξ − ′′ η 当 y(x)三次连续可微时, ) ( ) ( ) ( k k O hk y′′ ξ − y′′ η = 。忽略 ( ) 3 O hk 项，用 1 , k k+ y y 分 别近似 ( ), ( ) k k+1 y x y x ，得 [ ( , ) ( , )] 2 +1 = + k k + k+1 k+1 k k k f x y f x y h y y 称为梯形方法。 1.1.3 改进欧拉方法 在实际计算中，先用显式 Euler 方法所得的 k+1 y 为 (0) k+1 y ，再用梯形方法改进 一次 1 1 11 (, ) [ ( , ) ( , )], 2 0,1, , 1 k k kk k k kk k k y y hf x y h y y fx y fx y k n + + ++ ⎧ = + ⎪ ⎪ ⎨ =+ + ⎪ ⎪ = − ⎩ " 称为改进 Euler 方法，改进 Euler 方法还可写成： [ ( , ) ( , ( , ))] 2 1 k k k 1 k k k k k k k f x y f x y h f x y h y + = y + + + + 或 1 1 2 2 2 k h k h y y k k k + = k + + ) ( , 1 k k k = f x y ) ( , 2 1 1 k f x y h k = k+ k + k 1.1.4 单步法的局部截断误差、整体截断误差 设所用数值法为 ⎩ ⎨ ⎧ = + = − = + + ( , , , ), 0,1 , 1 1 1 0 y y h x y y h k n y k k kφ k k k k " α 设 y(x)是(1)的准确解,称 ) ( ) ( ) ( , ( ), ( ), k 1 k 1 k k k k k 1 hk y x y x h x y x y x + = + − − φ + τ 为此数值法在 k+1 x 的局部截断误差；称 ) ek = y(xk ) − yk (k = 0,1,",n 为此数值法在 k x 点的整体截断误差;如果对充分小的h 0, x a kh, > k = + 成立 max ( ) ( 1) 0 − ≤ ≥ ≤ ≤ y x y Ch p p k k k n 常数C 独立于h ,就称此方法是 p 阶方法