偏微分方程程序附件

合集下载

(整理)偏微分方程word电子讲义.

偏微分方程偏微分方程是一个非常广泛的课题，它包含分析的许多方面内容。

就我们现在的知识水平来说，我们只了解很少一点东西。

从十八世纪初开始，人们就开始结合物理、力学问题来研究偏微分方程，最早研究的几个方程是弦振动方程、热传导方程及调和方程，这部分理论已经被彻底地研究了，而且近乎完备，把它们称为偏微分方程的古典理论。

十八世纪后期在连续介质力学中研究流体的运动规律，在考虑流体的粘性时，描述运动规律的方程称为Navier-Stokes方程组，而在不计流体的粘性时，称为Euler方程组。

在此时期，描述弹性体运动规律的方程称为Saint Venant方程组。

到了十九、二十世纪，人们发现了描述电磁场运动规律的Maxwell方程组，描述微观粒子运动规律的Schrodinger方程及Dirac方程组，广义相对论中确定引力场的基本方程Einstein方程以及基本粒子规范场理论的基本方程Yang-Mills方程，在微分几何中研究极小曲面的极小曲面方程等等。

随着科学理论变得复杂，所提出的偏微分方程就愈多而且更加变化多端，可能出现的偏微分方程和方程组类型之多是出于想象的。

我们的目的是介绍现代偏微分方程理论中用到的一些技巧和方法。

众所周知，一本偏微分方程的书只能包括已有的基本材料的一小部分，因此我们必须作出选择，如何选择不是立足于逻辑基础上的，这种选择的主观性是相当明显的。

偏微分方程的内容是研究偏微分方程解的各种性质。

通常考虑以下问题1．对单个方程或方程组，应配以怎样的初值条件与边值条件使之具有解，这是解的存在性问题。

在研究解的存在性时，要明确解赖以存在的函数类。

2．解的唯一性或究竟有几个解，要明确使解为唯一的函数类。

3．解的正则性或光滑性。

是否为古典解、强解还是弱解？解具有几阶可微性？4．解的连续依赖性，必须明确是什么空间、什么范数实现的。

通常考虑的是解关于初、边值或关于方程系数，或在方程为线性时关于自由项的连续依赖性。

5．定解区域与影响区域。

(高等数学)偏微分方程

第十四章偏微分方程物理、力学、工程技术和其他自然科学经常提出大量的偏微分方程问题.由于实践的需要和一些数学学科（如泛函分析，计算技术）的发展，促进了偏微分方程理论的发展，使它形成一门内容十分丰富的数学学科.本章主要介绍一阶偏微分方程、线性方程组及二阶线性偏微分方程的理论.在二阶方程中，叙述了极值原理、能量积分及惟一性定理.阐明了一些解的性质和物理意义，介绍典型椭圆型、双曲型、抛物型方程的常用解法：分离变量法，基本解，格林方法,黎曼方法,势位方法及积分变换法.最后，扼要地介绍了有实用意义的数值解法：差分方法和变分方法.§1 偏微分方程的一般概念与定解问题[偏微分方程及其阶数] 一个包含未知函数的偏导数的等式称为偏微分方程.如果等式不止一个，就称为偏微分方程组.出现在方程或方程组中的最高阶偏导数的阶数称为方程或方程组的阶数.[方程的解与积分曲面] 设函数u 在区域D 内具有方程中所出现的各阶的连续偏导数，如果将u 代入方程后，能使它在区域D 内成为恒等式，就称u 为方程在区域D 中的解，或称正规解. ),,,(21n x x x u u = 在n +1维空间),,,,(21n x x x u 中是一曲面，称它为方程的积分曲面. [齐次线性偏微分方程与非齐次线性偏微分方程] 对于未知函数和它的各阶偏导数都是线性的方程称为线性偏微分方程.如()()()()y x f u y x c yuy x b x u y x a ,,,,=+∂∂+∂∂就是线性方程.在线性方程中，不含未知函数及其偏导数的项称为自由项，如上式的f (x,y ).若自由项不为零，称方程为非齐次的.若自由项为零，则称方程为齐次的.[拟线性方程与半线性方程] 如果一个方程，对于未知函数的最高阶偏导数是线性的，称它为拟线性方程.如()()()()()()0,,,,,,,,,,,,22222122211=+∂∂+∂∂+∂∂+∂∂∂+∂∂u y x c y uu y x b x u u y x a yu u y x a y x u u y x a x u u y x a就是拟线性方程，在拟线性方程中，由最高阶偏导数所组成的部分称为方程的主部.上面方程的主部为()()()22222122211,,,,,,yuu y x a y x u u y x a x u u y x a ∂∂+∂∂∂+∂∂如果方程的主部的各项系数不含未知函数，就称它为半线性方程.如()()()()0,,,,,,2222=∂∂+∂∂+∂∂+∂∂y yu y x d x y u y x c yu y x b x u y x a就是半线性方程.[非线性方程] 不是线性也不是拟线性的方程称为非线性方程.如1)()1(222=∂∂+∂∂+yux u u就是一阶非线性偏微分方程.[定解条件] 给定一个方程，一般只能描写某种运动的一般规律，还不能确定具体的运动状态，所以把这个方程称为泛定方程.如果附加一些条件（如已知开始运动的情况或在边界上受到外界的约束）后，就能完全确定具体运动状态，称这样的条件为定解条件.表示开始情况的附加条件称为初始条件，表示在边界上受到约束的条件称为边界条件.[定解问题] 给定了泛定方程（在区域D 内）和相应的定解条件的数学物理问题称为定解问题.根据不同定解条件，定解问题分为三类.1︒ 初值问题只有初始条件而没有边界条件的定解问题称为初值问题或柯西问题. 2︒ 边值问题只有边值条件而没有初始条件的定解问题称为边值问题.3︒ 混合问题既有边界条件也有初始条件的定解问题称为混合问题（有时也称为边值问题）.[定解问题的解] 设函数u 在区域D 内满足泛定方程，当点从区域D 内趋于给出初值的超平面或趋于给出边界条件的边界曲面时，定解条件中所要求的u 及它的导数的极限处处存在而且满足相应的定解条件，就称u 为定解问题的解.[解的稳定性] 如果定解条件的微小变化只引起定解问题的解在整个定义域中的微小变化，也就是解对定解条件存在着连续依赖关系，那末称定解问题的解是稳定的.[定解问题的适定性] 如果定解问题的解存在与惟一并且关于定解条件是稳定的，就说定解问题的提法是适定的.§2 一阶偏微分方程一、柯西-柯娃列夫斯卡娅定理[一阶偏微分方程的通解] 一阶偏微分方程的一般形式是0),,,,,,,,(2121=∂∂∂∂∂∂nn x ux u x u u x x x F或()0,,,,,,,211=n n p p p u x x F ，其中()n i x up ii ,,2,1 =∂∂=如解出p 1，可得：p 1 = f (x 1 , x 2 ,…, x n , u , p 2 ,…, p n )当方程的解包含某些“任意元素”(指函数)，如果适当选取“任意元素”时，可得方程的任意解(某些“奇异解”除外)，则称这样的解为通解.在偏微分方程的研究中，重点在于确定方程在一些附加条件(即定解条件)下的解，而不在于求通解.[一阶方程的柯西问题]()()⎪⎩⎪⎨⎧==∂∂=n x x n n x x u p p u x x x f x u,,|,,,,,,,22211011 ϕ 称为柯西问题，式中),,(2n x x ϕ为已知函数，对柯西问题有如下的存在惟一性定理.[柯西－柯娃列夫斯卡娅定理] 设 f ( x 1 , x 2 ,, x n , u , p 2 ,, p n ) 在点 ( x 10 , x 20 ,, x n 0 , u 0 , p 20 ,, p n 0 ) 的某一邻域内解析，而),,(2n x x ϕ在点( x 20 ,, x n 0 ) 的某邻域内解析，则柯西问题在点 ( x 10 ,, x n 0 ) 的某一邻域内存在着惟一的解析解.这个定理应用的局限性较大，因它要求f 及初始条件都是解析函数，一般的定解问题未必能满足这种条件.对高阶方程也有类似定理.二、一阶线性方程1．一阶齐次线性方程[特征方程∙特征曲线∙初积分(首次积分)] 给定一阶齐次线性方程在有些书中写作0),,,,,,,,,(121=∂∂∂∂∂∂nn x u x u t u u x x x t F()()0,,,,,,211211=∂∂++∂∂nn n n x u x x x a x u x x x a (1) 式中a i 为连续可微函数，在所考虑的区域内的每一点不同时为零(下同).方程组()n i ix x x a tx ,,,d d 21 = ( i = 1,2,, n ) 或()()()n n n n n x x x a x x x x a x x x x a x ,,,d ,,,d ,,,d 2121222111 === (2)称为一阶齐次线性偏微分方程的特征方程.如果曲线l : x i = x i (t ) ( i =1,2,, n )满足特征方程(2)，就称曲线l 为一阶齐次线性方程的特征曲线.如果函数ψ ( x 1 , x 2 ,, x n )在特征曲线),,2,1()(n i t x x i i ==上等于常数，即ψ ( x 1(t ) , x 2(t ) ,, x n (t ) ) = c就称函数ψ ( x 1, x 2,, x n )为特征方程(2)的初积分(首次积分). [齐次方程的通解]1o 连续可微函数u = ψ ( x 1, x 2,, x n ) 是齐次线性方程(1)的解的充分必要条件是： ψ ( x 1, x 2,, x n )是这个方程的特征方程的初积分.2o 设ψi ( x 1 , x 2 ,, x n ) ( i = 1,2,, n 1-) 是特征方程(2)在区域D 上连续可微而且相互独立的初积分(因此在D 内的每一点，矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂---n n n n n n x x x x x x x x x 121112221212111ψψψψψψψψψ 的秩为n 1-) ，则u = ω ( ψ1 ( x 1 , x 2 ,, x n ) ,, ψn -1 ( x 1 , x 2 ,, x n ) )是一阶齐次线性方程(1)的通解，其中ω为n 1-个变量的任意连续可微函数. [柯西问题] 考虑方程的柯西问题()()⎪⎩⎪⎨⎧==∂∂==∑n x x ni i n i x x u x u x x x a ,,|0,,,2121011 ϕ 式中ϕ ( x2 ,, x n )为已知的连续可微函数.设ψi ( x 1 , x 2 ,, x n ) ( i = 1,2,, n 1-) 为特征方程的任意n 1-个相互独立的初积分，引入参变量 i ψ (1,,2,1-=n i )，从方程组()()()⎪⎪⎩⎪⎪⎨⎧===--120112201212011,,,,,,,,,n n n n n x x x x x x x x x ψψψψψψ 解出x 2 ,, x n 得()()⎪⎩⎪⎨⎧==--12112122,,,,,,n n nn x x ψψψωψψψω 则柯西问题的解为u = ϕ ( ω2 ( ψ1 , ψ2 ,, ψn -1 ) ,, ωn ( ψ1 , ψ2 ,, ψn -1 ) )2．非齐次线性方程它的求解方法与拟线性方程相同.三、一阶拟线性方程一阶拟线性方程为()()∑==∂∂ni n i n i u x x x R x uu x x x a 12121,,,,,,,, 其中a i 及R 为x 1 , x 2 ,, x n , u 的连续可微函数且不同时为零. [一阶拟线性方程的求解和它的特征方程]()()⎪⎩⎪⎨⎧===u x x x R t un i u x x x a t x n n i i,,,,d d ),,2,1(,,,,d d 2121 或()()()u x x R uu x x a x u x x a x n n n n n ,,,d ,,,d ,,,d 11111 ===为原拟线性方程的特征方程.如果曲线l : x i = x i (t ) ( i =1,2,, n ) , u = u (t ) 满足特征方程，则称它为拟线性方程的特征曲线.设 ψi ( x 1 ,, x n ,u ) ( i = 1,2,, n ) 为特征方程的n 个相互独立的初积分，那末对于任何连续可微函数ω，ω ( ψ1 ( x 1,, x n , u ) , ψ2 ( x 1,, x n , u ) ,, ψn ( x 1,, x n , u ) ) = 0都是拟线性方程的隐式解.[柯西问题] 考虑方程的柯西问题()()()⎪⎩⎪⎨⎧==∂∂==∑n x x ni n i ni x x u u x x x R x u u x x x a ,,|,,,,,,,,212121011 ϕ ϕ为已知的连续可微函数.设 ψ1 ( x 1 , x 2 ,, x n , u ) ,, ψn ( x 1 , x 2 ,, x n , u ) 为特征方程的n 个相互独立的初积分，引入参变量 n ψψψ,,,21 , 从()()()⎪⎪⎩⎪⎪⎨⎧===nn n n n u x x x u x x x u x x x ψψψψψψ,,,,,,,,,,,,2012201212011解出 x 2 ,, x n , u()()()⎪⎪⎩⎪⎪⎨⎧===n n n n n u x x ψψψωψψψωψψψω,,,,,,,,,21212122 则由()()()()()()()0,,,,,,,,,,,,,,,,,,,,,,2121221221121=-≡n n n n n n u x x x u x x x u x x x V ψψψωψψψωϕψψω给出柯西问题的隐式解.四、一阶非线性方程[完全解·通解·奇异解] 一阶非线性方程的一般形式为()()n i x up p p p u x x x F ii n n ,,2,10,,,,,,,,2121 =∂∂== 若一阶偏微分方程的解包含任意n 个独立的常数，则称这样的解为完全解(全积分). 若V ( x 1, x 2 ,, x n , u , c 1 , c 2,, c n ) = 0为方程的完全解，从()n i c VV i,,2,10,0 ==∂∂= 消去c i ，若得一个解，则称它为方程的奇异解(奇积分).以两个独立变量为例说明完全解与通解、奇异解的关系，设方程()yzq x z p q p z y x F ∂∂=∂∂==,,0,,,,有完全解V (x ,y ,z ,a ,b )=0 ( a ,b 为任意常数),则方程等价于从方程组()⎪⎪⎩⎪⎪⎨⎧=∂∂+∂∂=∂∂+∂∂=0,00,,,,q z Vy V p z V x V b a z y x V 消去a ,b 所得的方程.利用常数变易法把a ,b 看作x , y 的函数，将V (x ,y ,z ,a ,b )=0求关于x , y 的偏导数，得00=∂∂⋅∂∂+∂∂⋅∂∂+∂∂+∂∂=∂∂⋅∂∂+∂∂⋅∂∂+∂∂+∂∂ybb V y a a V q z V y V xbb V x a a V p z V x V那末0,0=∂∂⋅∂∂+∂∂⋅∂∂=∂∂⋅∂∂+∂∂⋅∂∂yb b V y a a V x b b V x a a V 与V=0联立可确定a ,b .有三种情况：1︒ 0≡∂∂≡∂∂bVa V ，将其与V (x ,y ,z ,a ,b )=0联立可确定不含任意常数的奇异解. 2︒ 如0=∂∂=∂∂=∂∂=∂∂yb x b y a x a ，即回到完全解. 3︒ 当0/,0/≡∂∂≡∂∂b Va V 时，必有()()0,,=∂∂y x b a ，这时，如果不属于情形2︒ ，则a 与b 存在函数关系：b=ω(a )，这里ω为任意可微函数，并从方程V (x ,y ,z ,a ,b )=0和()∂∂∂∂ωV a Vba +'=0消去a ,b ，可确定方程的通解.定理偏微分方程的任何解包含在完全解内或通解内或奇异解内. [特征方程·特征带·特征曲线·初积分] 在一阶非线性方程：()F x x x u p p p n n 12120,,,,,,,, =中，设F 对所有变量的二阶偏导数存在且连续，称()n i uFp x F t p p F p t u p Ft x i i i ni iii i ,,2,1)(d d d d ,1 =∂∂+∂∂-=∂∂=∂∂=∂∂∑=或u F p x F p u F p x F p p Fp up F x p F xp F x n nnni i i nn ∂∂+∂∂-==∂∂+∂∂-=∂∂=∂∂==∂∂=∂∂∑=d d d d d d 11112211为非线性方程的特征方程.设特征方程的解为x i =x i (t ), u=u (t ), p i =p i (t ) (i =1,2,…,n )称它为非线性方程的特征带.在x 1,x 2,, x n ,u 空间的曲线x i =x i (t ), u=u (t ) (i=1,2,…,n )称为非线性方程的特征曲线.如果函数()n n p p p u x x x G ,,,,,,,,2121 在特征方程的任一解x i =x i (t ) (i =1,2,, n ), u=u (t ), p i =p i (t ) (i =1,2,, n )上等于常数，即()()()()()()()()G x t x t x t u t p t p t p t C n n 1212,,,,,,,, =那末函数()n n p p p u x x x G ,,,,,,,,2121 称为特征方程的初积分.[求完全解的拉格朗日－恰比方法] 考虑两个变量的情况.对于方程F (x ,y ,z ,p ,q )=0，选择使雅可比式()()0,,≠∂∂q p G F 的一个初积分G (x ,y ,z ,p ,q ).解方程组()()F x y z p q G x y z p q a,,,,,,,,==⎧⎨⎪⎩⎪0(a 为任意常数) 得p (x ,y ,z ,a )及q (x ,y ,z ,a ).则方程d z=p d x+q d y的通解V (x ,y ,z ,a ,b )=0(b 是积分d z=p d x+q d y 出现的任意常数)就是方程F (x ,y ,z ,p ,q )=0的完全解.例求方程()z p q x y 22222+=+的完全解.解方程的特征方程为()()()qy x z y qp q p z x p q p z z q z y p z x 22222222222d 22d 2d 2d 2d +-=+-=+== 这里成立zpxx p z z p d d d =+ 所以特征方程的一个初积分为z 2p 2 －x 2 .解方程组 ()()z p q x y z p x a22222222+-+=-=⎧⎨⎪⎩⎪ (a 为任意常数) 得 p a x zq y az=+=-22, 积分微分方程dz a x zdx y azdy =++-22 得完全解z x x a y y a a x x a y y ab 22222=++-++++-+ln(b 为任意常数)[某些容易求完全解的方程] 1︒ 仅含p ,q 的方程F (p ,q )=0G =p 是特征方程的一个初积分.从F (p ,q )=0与p=a (a 为任意常数)得q=ψ(a )，积分d z=a d x+ψ(a )d y得完全解z=ax+ψ(a )y+b (b 为任意常数)2︒ 不显含x ,y 的方程F (z ,p ,q )=0 特征方程为zFqqz F p p q F q p F p z q F y p F x ∂∂-=∂∂-=∂∂+∂∂=∂∂=∂∂d d d d d 因此q d p-p d q =0，显然G qp=为一个初积分，由F (z ,p ,q )=0，q=pa (a 为任意常数)解得p=ψ(z ,a ).于是由d z=ψ(z ,a )d x+a ψ(z ,a )d y得()⎰++=b ay x a z z,d ψ (b 为任意常数)可确定完全解.3︒ 变量分离形式的方程()f x p i i i i n,=∑=10特征方程为n n n n i i iin n n x f p x f p p f p z p f x p f x ∂∂-==∂∂-=∂∂=∂∂==∂∂∑=d d d d d 1111111 可取初积分G i =f i (x i ,p i ) , (i =1,2,, n ).从f i (x i ,p i )=a i (i =1,2,, n )解出p i =ϕi (x i ,a i )得完全解()∑⎰=+=ni i i i i b x a x z 1d ,ϕ式中a i ,b 为任意常数，且a i i n=∑=10.[克莱罗方程] 方程()z p x f p p p i i n i n=+=∑121,,,称为克莱罗方程，其完全解为()z c x f c c c i i n i n=+=∑121,,,对c i 微分得x fc i i=-∂∂ (i =1,2,…,n ) 与完全解的表达式联立消去c i 即得奇异解.例求方程z －xp －yq －pq =0的完全解和奇异解. 解这是克莱罗方程，它的完全解是z=ax+by+ab对a,b 微分，得x=－b,y=－a ，消去a ,b 得奇异解z=－xy[发甫方程] 方程P (x,y,z )d x+Q (x,y,z )d y+R (x,y,z )d z=0 (1)称为发甫方程，如果P,Q,R 二次连续可微并满足适当条件，那末方程可积分.如果可积分成一关系式时，则称它为完全可积.1︒ 方程完全可积的充分必要条件当且仅当P,Q,R 满足条件0)()()(=∂∂-∂∂+∂∂-∂∂+∂∂-∂∂yP x Q R x R z P Q z Q y R P (2) 时，存在一个积分因子μ(x,y,z )，使d U 1=μ(P d x+Q d y+R d z )从而方程的通解为U 1(x,y,z )=c特别，当0,0,0=∂∂-∂∂=∂∂-∂∂=∂∂-∂∂yP x Q x R z P z Q y R 时，存在一个函数U (x,y,z )满足 zU R y U Q x U P ∂∂=∂∂=∂∂=,,从而 d U=P d x+Q d y+R d z 所以方程的通解为U (x,y,z )=c所以完全可积的发甫方程的通解是一单参数的曲面族.定理设对于发甫方程(1)在某区域D 上的完全可积条件(2)成立，则对D 内任一点M (x,y,z )一定有方程的积分曲面通过，而且只有一个这样的积分曲面通过. 2︒ 方程积分曲面的求法设完全可积条件(2)成立.为了构造积分曲面，把z 看成x,y 的函数(设R (x,y,z )≠0)，于是原方程化为y RQ x R P z d d d --=由此得方程组()()()()⎪⎪⎩⎪⎪⎨⎧≡-=∂∂≡-=∂∂4,,3,,11z y x Q R Q y z z y x P R P xz发甫方程(1)与此方程组等价.把方程(3)中的y 看成参变量，积分后得一个含有常数 c 的通解 ()cy x z ~;,ϕ= 然后用未知函数()~cy 代替常数 c ，将()()z x y c y =ϕ,;~代入方程(4)，在完全可积的条件下，可得()~cy 的一个常微分方程，其通解为 ()()~,cy y c =ψ c 为任意常数，代回()()z x y cy =ϕ,;~中即得发甫方程的积分曲面 z=ϕ(x,y,ψ(y,c ))由于发甫方程关于x,y,z 的对称性，在上面的讨论中，也可把x 或y 看成未知函数，得到同样的结果.例求方程yz d x+2xz d y+xy d z=0的积分曲面族.解容易验证完全可积条件成立，显然存在一个积分因子μ=1xyz，用它乘原方程得 0d d 2d =++zz y y x x 积分后得积分曲面族xy 2z=c也可把方程化为等价的方程组⎪⎪⎩⎪⎪⎨⎧-=∂∂-=∂∂y z yz x z xz 2 把y 看成参变量，积分xzx z -=∂∂得通解 zx c= 用未知函数()~cy 代替 c ，将()y c zx ~=代入方程y z y z 2-=∂∂得 ()()yy cy y c ~2d ~d -= 积分后有()~cy c y =2所以原方程的积分曲面族是xy 2z=c五、一阶线性微分方程组[一阶线性偏微分方程组的一般形式] 两个自变量的一阶线性方程组的形式是()n i F u C x u B t u A i n j j ij n j n j jij j ij ,,2,10111 ==++∂∂+∂∂∑∑∑=== 或()n i f u b x u a t u i n j j ij n j j ij i,,2,1011 ==++∂∂+∂∂∑∑== (1) 其中A ij ,B ij ,C ij ,F i ,a ij ,b ij ,f i 是(x,t )的充分光滑函数. [特征方程·特征方向·特征曲线]⎩⎨⎧=≠==-j i j i t xa ij ij ij ,1,0,0)d d det(δδ称为方程组(1)的特征方程.在点(x,t )满足特征方程的方向txd d 称为该点的特征方向.如果一条曲线l ，它上面的每一点的切线方向都和这点的特征方向一致，那末称曲线l 为特征曲线. [狭义双曲型方程与椭圆型方程] 如果区域D 内的每一点都存在n 个不同的实的特征方向，那末称方程组在D 内为狭义双曲型的.如果区域D 内的每一点没有一个实的特征方向，那末称方程组在D 内为椭圆型的. [狭义双曲型方程组的柯西问题] 1︒ 化方程组为标准形式——对角型因为det(a ij -δij λ)=0有n 个不同的实根λ1(x,t ) ,, λn (x,t )，不妨设),(),(),(21t x t x t x n λλλ<<<那末常微分方程()()n i t x txi ,,2,1,d d ==λ 的积分曲线l i (i =1,2,…,n )就是方程组(1)的特征曲线. 方程()()aijk ij k i i n-==∑λδλ1的非零解(λk (1) ,, λk (n ))称为对应于特征方向λk 的特征矢量. 作变换()()n i u v nj jj i i ,,2,11==∑=λ可将方程组化为标准形式——对角型()()()()n i t x v t x a x v t x t v i nj j ij ii i ,,2,1,,,1=+=∂∂+∂∂∑=βλ 所以狭义双曲型方程组可化为对角型，而一般的线性微分方程组(1)如在区域D 内通过未知函数的实系数可逆线性变换可化为对角型的话，(此时不一定要求 λi 都不相同)，就称这样的微分方程组在D 内为双曲型的. 2︒ 对角型方程组的柯西问题考虑对角型方程组的柯西问题()()()()()()n i x x v t x v t x a x v t x tv i inj i j ij i i i,,2,10,,,,1 =⎪⎩⎪⎨⎧=+=∂∂+∂∂∑=ϕβλ ϕi (x )是[a,b ]上的连续可微函数.设αij ,βi ,λi 在区域D 内连续可微，在D 内可得相应的积分方程组()()()n i tv x t x v il i n j j ij i i i ,,2,1d ,~1 =⎥⎦⎤⎢⎣⎡++=⎰∑=βαϕ 式中 l i 为第i 条特征曲线l i 上点(x,t )与点(x i ,0)之间的一段，(x i ,0)为l i与x 轴上[a,b ]的交点.上式可以更确切地写为()()[]()[]()[]()[]⎰∑⎭⎬⎫⎩⎨⎧+⋅+==t n j i i i j i ij i i i t x x t x x v t x x a t x x t x v 01d ,,,,,,,,,0,,,τττβττττϕ(i =1,2,, n )式中x i =x i (x ︒,t ︒,t )为过点(x ︒,t ︒)的第i 条特征曲线，利用逐次逼近法可解此积分方程.为此令()()()[]()()()()[]()[]()()[]()[]()()()()[]()[]()()[]()[]()n i t x x t x x v t x x a t x x t x v n i t x x t x x v t x x a t x x t x v n i t x x t x v i i tnj i k j i ij i i k ii i tnj i j i ij i i ii i i ,,2,1d ,,,,,,,,,0,,,,,2,1d ,,,,,,,,,0,,,,,2,10,,,}{}{01101010=+⋅+==+⋅+===⎰∑⎰∑=-=τττβττττϕτττβττττϕϕ序列{v i (k )} (k =0,1,2 ,)一致收敛于积分方程的连续可微解v i (x,t ) (i =1,2,, n )，这个v i (x,t )也就是对角型方程组的柯西问题的解.设在区域D 内对角型方程组的柯西问题的解存在，那末解与初值有下面的关系：(i) 依赖区间：过D 中任意点M (x,t )作特征曲线l 1,l n ，交x 轴于B,A ，称区间[A,B ]为M 点的依赖区间(图14.1(a ))，解在M 点的值由区间[A,B ]的初值确定而与[A,B ]外的初值无关. (ii) 决定区域：过点A,B 分别作特征曲线l n ,l 1,称l n ,l 1 与区间[A,B ]围成的区域D 1为区间[A,B ]的决定区域(图14.1(b ))，在区域D 1中解的值完全由[A,B ]上的初值决定.(iii) 影响区域：过点A,B 分别作特征曲线l 1,l n ，称l 1,l n 与[A,B ]围成的区域D 2为区间[A,B ]的影响区域(图14.1(c )).特别当区间[A,B ]缩为一点A 时，A 点的影响区域为D 3(图14.1(d )).在区域D 2中解的值受[A,B ]上的初值影响，而在区域D 2外的解的值则不受[A,B ]上的初值影响.图14.1[线性双曲型方程组的边值问题] 以下列线性方程组来说明：()⎪⎪⎩⎪⎪⎨⎧<++=∂∂+∂∂++=∂∂+∂∂2122221111λλλλc v b u a x v t v c v b u a xu t u (1) 1︒ 第一边值问题(广义柯西问题) 设在平面(x,t )上给定曲线段⋂AB ，它处处不与特征方向相切.过A,B 分别引最左和最右的特征曲线l 1及l 2.要求函数u (x,t ),v (x,t )在⋂AB ，l 1及l 2围成的闭区域D 上满足方程组，且在⋂AB 上取给定的函数值(图14.2(a )).2︒ 第二边值问题(古沙问题) 设l 1是过P 点的第一族特征线，l 2是第二族特征线，在l 1的一段PA 上给定v (x,t )的数值，在l 2的一段PB 上给定u (x,t )的数值，过A 点作第二族特征线，过B 点作第一族特征线相交于Q .求在闭区域PAQB 上方程组的解(图14.2(b )).3︒ 第三边值问题设AB 为非特征曲线的曲线弧，AC 为一特征线弧，且在AB 与AC 之间不存在过A 点的另外特征曲线,过C 点作第二族特征线与过B 点的第一族特征线交于E 点，在AC 上给定v (x,t )的数值，在AB 上给定u (x,t )的数值，求ACEBA 所围成的闭区域D 上的方程组的解(图14.2(c)).图14.2[边值问题的近似解——特征线法] 以上定解问题，可用逐步逼近法求解，也可用特征线法求解的近似值.以第一边值问题为例说明.在曲线AB 上取n 个分点A 1,A 2,, A n ，并记A 为A 0，B 为A n +1，过A 0按A 0的第二特征方向作直线与过A 1按A 1的第一特征方向作直线相交于B 0；过A 1按A 1第二特征方向作直线与过A 2按A 2的第一特征方向作直线相交于B 1 ,最后得到B n (图14.3).用如下的近似公式来确定方程组(1)的解u (x,t ),v (x,t )在B i (i =0,1,2,…,n )的数值：()()()()()()(){}()[]()()()()()()(){}()[]u B u A B A a A u A b A v A c A A v B v A B A a A u A b A v A c A A i i i i i i i i i i i i i i i i i i i i -=++⨯+-=++⨯+⎧⎨⎪⎩⎪+++++++--11111111112122212121211λλ图14.3于是在一个三角形网格的节点上得到u,v 的数值.再经过适当的插值，当n 相当大，A i 、A i +1的距离相当小时，就得到所提问题的足够近似的解.[特殊形式的拟线性方程组——可化约系统] 一般的拟线性方程组的问题比较复杂，目前研究的结果不多，下面介绍一类特殊形式的拟线性方程组——可化约系统.如果方程组⎪⎪⎩⎪⎪⎨⎧=∂∂+∂∂+∂∂+∂∂=∂∂+∂∂+∂∂+∂∂0022221111x v D t v C x u B tu A xv D t v C x u B t uA 中所有的系数只是u,v 的函数，称它为可化约系统. 考虑满足条件()()0,,≠∂∂t x v u 的方程组的解u=u (x,t ),v=v (x,t ).x,t 可以表示成u,v 的函数，且()()()()()()()()v u t x u t x v v u t x u x t v v u t x v tx u v u t x v x t u ,,,,,,,,,,∂∂∂∂=∂∂∂∂∂∂-=∂∂∂∂∂∂-=∂∂∂∂∂∂=∂∂ 原方程化为⎪⎪⎩⎪⎪⎨⎧=∂∂+∂∂-∂∂-∂∂=∂∂+∂∂-∂∂-∂∂0022221111u t D u x C v t B vx A ut D u x C v t B v xA 这是关于自变量u,v 的线性方程组.这样就把求拟线性方程组满足()()0,,≠∂∂t x v u 的解，化为解线性方程组的问题.而此线性方程组满足条件()()0,,≠∂∂v u t x 的解，在(x,t )平面上的象即为原来拟线性方程组的解.§3 二阶偏微分方程一、二阶偏微分方程的分类、标准形式与特征方程考虑二阶偏微分方程()0),,,,,,(111,2=∂∂∂∂+∂∂∂∑=nnnj i j i ij x u x u u x x F y x u x a (1) 式中a ij (x )=a ij (x 1,x 2,…,x n )为x 1,x 2,…,x n 的已知函数.[特征方程·特征方向·特征曲面·特征平面·特征锥面]代数方程()01,=∑=nj i jiijaa x a称为二阶方程(1)的特征方程；这里a 1,a 2,…,a n 是某些参数，且有012≠∑=ni i a .如果点x ︒=(x 1︒,x 2︒,…,x n ︒)满足特征方程，即()01,o =∑=nj i jiijaa x a则过x ︒的平面()01o=-∑=nk kk k x x a 的法线方向l :(a 1,a 2,…,a n )称为二阶方程的特征方向；如果一个(n 1-)维曲面，其每点的法线方向都是特征方向，则称此曲面为特征曲面；过一点的(n 1-)维平面，如其法线方向为特征方向，则称这个平面为特征平面，在一点由特征平面的包络组成的锥面称为特征锥面.[n 个自变量方程的分类与标准形式] 在点P (x 1︒,x 2︒,…,x n ︒)，根据二次型()∑=nj i jinijaa x x x a 1,o o 2o 1,,, (a i 为参量)的特征根的符号，可将方程分为四类：(i) 特征根同号，都不为零，称方程在点P 为椭圆型.(ii) 特征根都不为零，有n 1-个具有同一种符号，余下一个符号相反，称方程在点P 为双曲型.(iii) 特征根都不为零，有m n -个具有同一种符号(n >m >1)，其余m 个具有另一种符号，称方程在点P 为超双曲型.(iv) 特征根至少有一个是零，称方程在点P 为抛物型.若在区域D 内每一点方程为椭圆型，双曲型或抛物型，则分别称方程在区域D 内是椭圆型、双曲型或抛物型.在点P 作自变量的线性变换可将方程化为标准形式：椭圆型：∑==+∂∂ni ix u1220Φ双曲型：∑==+∂∂-∂∂n i ix ux u 22120Φ超双曲型：()10112222>>=+∂∂-∂∂∑∑=+=m n x ux u m i nm i ii Φ抛物型：()00122>=+∂∂∑-=m x umn i iΦ 式中Φ为不包含二阶导数的项.[两个自变量方程的分类与标准形式] 方程的一般形式为0,,,,222222122211=⎪⎪⎭⎫⎝⎛∂∂∂∂+∂∂+∂∂∂+∂∂y u x u u y x F y u a y x u a x u a (2) a 11,a 12,a 22为x ,y 的二次连续可微函数，不同时为零. 方程a 11d y 22-a 12d x d y +a 22d x 2=0称为方程(2)的特征方程.特征方程的积分曲线称为二阶方程(2)的特征曲线. 在某点P (x 0,y 0)的邻域D 内，根据Δ=a 122－a 11a 12的符号将方程分类：当Δ>0时，方程为双曲型；当Δ=0时，方程为抛物型；当Δ<0时，方程为椭圆型.在点P 的邻域D 内作变量替换，可将方程化为标准形式：（i ）双曲型：因Δ>0，存在两族实特征曲线11),(c y x =ϕ，22),(c y x =ϕ，作变换),(1y x ϕξ=，),(2y x ϕη=和,,ηηξ-=+=s t s 方程化为标准形式),,,,(2222tus u u t s t u s u ∂∂∂∂=∂∂-∂∂Φ或),,,,(12ηξηξΦηξ∂∂∂∂=∂∂∂uu u u （ii ）抛物型: 因Δ=0，只存在一族实的特征曲线c y x =),(ϕ，取二次连续可微函数),(y x ψ，使0),(),(≠∂∂y x ψϕ，作变换),(y x ϕξ=，),(y x ψη=，方程化为标准形式),,,,(222ηξηξΦη∂∂∂∂=∂∂uu u u （iii ）椭圆型：因Δ<0，不存在实特征曲线，设c y x i y x y x =+=),(),(),(21ϕϕϕ为11221121212d d a a a a a x y -+=的积分，y x ϕϕ,不同时为零，作变量替换),(1y x ϕξ=，),(2y x ϕη=,方程化为标准形式),,,,(32222ηξηξΦηξ∂∂∂∂=∂∂+∂∂uu u u u二、极值原理·能量积分·定解问题的惟一性定理椭圆型方程、抛物型方程的极值原理及双曲型方程的能量守恒原理是相应方程的解所具有的最基本性质之一，在定解问题的研究中起着重要的作用. [椭圆型方程的极值原理与解的惟一性定理]1︒ 极值原理设D 为n 维欧氏空间E n 的有界区域，S 是D 的边界，在D 内考虑椭圆型方程()()()()x x x x f u c x ub x x u a Lu ni i i n j i j i ij =+∂∂+∂∂∂≡∑∑==11,2式中a ij (x ),b i (x ),c (x ),f (x )在D 上连续，c (x )≤0且二次型()∑=nj i j i ij a a a 1,x 正定，即存在常数μ>0，对任意x D ∈和任意的a i 有()∑∑==≥ni i nj i jiija aa a 121,μx定理1 设u (x )为D 内椭圆型方程的解，它在D 内二次连续可微，在D 上连续，且不是常数，如f (x )≤0（或f (x )≥0），则u (x )不能在D 的内点取非正最小值（或非负最大值）. 如果过边界S 上的任一点P 都可作一球，使它在P 点与S 相切且完全包含在区域D 内，则有定理2 设u (x )为椭圆型方程在D 内二次连续可微，在D 上连续可微的解，且不是常数，并设f (x )≤0（或f (x )≥0）.若u (x )在边界S 上某点M 处取非正最小值（或非负最大值），只要外法向导数错误！未定义书签。

c语言偏微分方程函数包,[转载]一个求解偏微分方程问题的C程序

c语⾔偏微分⽅程函数包,[转载]⼀个求解偏微分⽅程问题的C程序#include#include#define vareps 0.02#define gamma 0.1#define pi 3.1415926535898#define eps 1.0e-12#define x1 0#define x2 2*pi//参数取值double norm(double **a,int N,double h){double b=0;int i,j;for(i=0;ifor(j=0;jb+=a[i][j]*a[i][j]*h*h;b=sqrt(b);return b;}//误差计算，其中空间⽹格数N，步长hvoid main(){double **u0,**u,**uh,**e,*x,*y,r,h,tau,T;int i,j,k,K,N;FILE *fp=fopen("data.m","w");FILE *fp0=fopen("data1.m","w");printf("Printf T:n");scanf("%lf",&T);//输⼊时间参数Tprintf("Printf the length of tau:n");scanf("%lf",&tau);//时间步长tauprintf("Printf N:n");scanf("%d",&N);//空间⽹格数NK=(int)(T/tau)+1;h=(x2-x1)/N;//空间步长hN++;r=tau/(h*h);//⽹⽐ruh=new double*[N];for(i=0;iuh[i]=new double[N];u0=new double*[N];for(i=0;iu0[i]=new double[N];u=new double*[N];for(i=0;iu[i]=new double[N];e=new double*[N];for(i=0;ie[i]=new double[N];x=new double[N];y=new double[N];for(i=0;i{x[i]=x1+i*h;y[i]=x[i];}for(i=0;ifor(j=0;jif((x[i]-pi)*(x[i]-pi)+(y[j]-pi)*(y[j]-pi)<0.25) u[i][j]=1;else u[i][j]=-1;fprintf(fp0,"x=[");for(i=0;ifprintf(fp0,"%lft",x[i]);fprintf(fp0,"];ny=[");for(i=0;ifprintf(fp0,"%lft",y[i]);fprintf(fp0,"];nz=[");for(i=0;i{for(j=0;jfprintf(fp0,"%et",u[i][j]);fprintf(fp0,"n");}fprintf(fp0,"];nfigure,surf(x,y,z)nfigure,contour(x,y,z)n");//输出初始条件for(k=1;k{for(i=0;ifor(j=0;juh[i][j]=u[i][j];while(1){for(i=0;ifor(j=0;ju0[i][j]=u[i][j];u[0][0]=(uh[0][0]/tau/gamma+2/h/h*(u[0][1]+u[1][0])-u[0][0]*u[0][0]*u[0][0]/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps);u[N-1][0]=(uh[N-1][0]/tau/gamma+2/h/h*(u[N-1][1]+u[N-2][0])-u[N-1][0]*u[N-1][0]*u[N-1][0]/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps);u[0][N-1]=(uh[0][N-1]/tau/gamma+2/h/h*(u[0][N-2]+u[1][N-1])-u[0][N-1]*u[0][N-1]*u[0][N-1]/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps);u[N-1][N-1]=(uh[N-1][N-1]/tau/gamma+2/h/h*(u[N-2][N-1]+u[N-1] [N-2])-u[N-1][N-1]*u[N-1][N-1]*u[N-1][N-1]/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps);for(i=1;i{ u[i][0]=(uh[i][0]/tau/gamma+1/h/h*(2*u[i][1]+u[i-1][0]+u[i+1][0]) -u[i][0]*u[i][0]*u[i][0]/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps);u[i][N-1]=(uh[i][N-1]/tau/gamma+1/h/h*(2*u[i][N-2]+u[i-1][N-1]+u[i+1][N-1])-u[i][N-1]*u[i][N-1]*u[i][N-1]/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps);u[0][i]=(uh[0][i]/tau/gamma+1/h/h*(u[0][i+1]+u[0][i-1]+2*u[1][i])-u[0][i]*u[0][i]*u[0][i]/vareps/vareps)/(1/tau/gamma+4/h/h-1/vareps/vareps);u[N-1][i]=(uh[N-1][i]/tau/gamma+1/h/h*(u[N-1][i+1]+u[N-1][i-1]+2*。

(完整版)偏微分方程课程教学大纲

高维波动方程Cauchy问题
4
上课
作业
掌握思想
测验一
波动方程的有限传播速度
2
上课
作业
掌握波动方程的特性
一维波动方程初边值问题
4
上课
作业
掌握分离变量法
波动方程的能量方法
4
上课
作业
掌握PDE能量方法
期中考试
热方程的Cauchy问题
4
上课
作业
掌握Fourier变换及基本解思想
热ቤተ መጻሕፍቲ ባይዱ程的初边值问题
2
上课
作业
掌握分离变量法
在该课程中，我们将重点讲解刻画粒子输运、波的传播、温度变化及物理场论的输运方程、波动方程、热方程及Laplace方程的理论与方法，从而了解双曲型、抛物型及椭圆型三类偏微分方程的特性及基本分析思路,为进一步学习与研究打下一定的基础。
＊课程简介(Description)
In this course，we shall teach students how to formulate the partial differential equation models based on certain elementary laws in physics and mechanics, then purpose some analytic methods for studying qualitative and quantitative properties of solutions to several important kinds of partial differential equations,which could inspire students to understand some basic idea of modern methods and theories of partial differential equations。 By solving the explicit solutions of PDEs， it can explain some important phenomena in physics and mechanics, which could inspire the students to study further。

第一章偏微分方程概论

u0 , 0是由初始条件确定的2个常数。
1.1.2 一些典型的常微分方程
一、可分离变量的方程
具有如下形式：
d y f (x) g(y) 可转化为 d x
1 dy f (x) g( y) dx
两边对x积分（如果可能的话）
1 d y d x f( xd )x g ( y ) d x
xh
0
偏微分方程的导出与定解
且通解为
u = f (x - at) + g (x + at)
其中f与g是任意两个具有连续二阶导数的函数。并由初始条件，就得到下面弦振动的达朗贝尔（D′Alembert）公式
x a t u ( x a t )u ( x a t )1 0 0 u ud ( ) 1 2 x a t 2 2 a u
2 2 2 2 u 2 u u u a 2 2 2 2 t x y z
偏微分方程的导出与定解
1.2.3 初边值问题对于最典型的求解问题是初始值问题——柯西问题即：求波动方程的解 u ，使其满足初始条件
u(0 , x, y, z) u , y, z) 0(x u (0 , x, y, z) u , y, z) 1(x t
并由初始条件并由初始条件就得到下面弦振动的达朗就得到下面弦振动的达朗偏微分方程的导出与定解贝尔dalembert公式001u??1?2?2xat?xat??uxatuxatuud?a?????20???高维情形把xyz记xx1x2x3??1?2?3利用傅立叶变换fourier偏微分方程的导出与定解其中?x?1x1?2x2?3x3123123?f?fxxxeix?dxdxdx?????????????????且当f满足一定条件时有fourier逆变换12313331?ffx???2ddd????????????????????偏微分方程的导出与定解另外有123123?f?ix?iifi?xxxedxdxdxx?????????????????????222222222uuuuatxyz?????????????????对于下面方程利用fourier变换00xyzxyzuu????偏微分方程的导出与定解??222222312??duaudt???????10xyzxyzuut?0100??????ttduuuudt????变成解常微分方程的初值问题解得1230123??????cos?sinsinutuataatt??????偏微分方程的导出与定解1123??ua???其中做fourier逆变换得泊松poisson公式222231???????01111101221111utx4?4???114?4?ttatatttatuxatlds?atuxatlds?taaatuxldsatuxldstattat?????????????????????????????????????????偏微分方程的导出与定解其中ds1dsat是球面l1lat的面积元素

偏微分方程解的几道算例(差分、有限元)-含matlab程序(

《偏微分方程数值解》上机报告实验内容1:分别用向前差分格式、向后差分格式及六点对称格式，求解下列问题：222, 01, 0, (0, (1, 0,1, (, 0 sin( (1.u u x t t x u t u t t u x x x x ? ??=+<<>? ?n?? ==>?? =+- ? x 方向0.1h =, t 方向0.01 T在0.25t =时观察数值解与精确解2sin( (1 u e x x x 的误差+.-(—算法描述：匚实验结果：1.误差的数值解结果数值对比(A向前差分格式”程序：>>forward(0.1,0.01,0.25Current plot heldans =0.00000.00270.00510.00700.00820.00870.00820.00700.00510.00270.0000(B向后差分格式”程序：>>back(0.1,0.01,0.25Curre nt plot heldans =0.0000-0.0037-0.0071-0.0097-0.0114-0.0120-0.0114-0.0097-0.0071 -0.00370.0000(C六点差分格式”程序：>>six(0.1,0.01,0.25Current plot heldans =0.0000-0.0005-0.0009-0.0013-0.0015-0.0016-0.0015-0.0013-0.0009-0.00050.0000注:这里的"误差"=精确解-数值解.2.精确解与数值解结果图像对比向前差分格式注:曲线表示精确解，"o"表示数值解（t=0.25时.六点差分格式0^60 250150<K注:曲线表示精确解,"o"表示数值解（t=0.25时. “后差分格式a 7j □? 0.30^6注:曲线表示精确解，"0"表示数值解（t=0.25时.（三结果分析通过（一,（二，我们检验了三种方法都能很好的求解此一维热传导方程，其中明显能发现六点对称格式”的误差更小。

偏微分方程数值解流程

偏微分方程数值解流程1.网格划分：将求解域划分为网格，这是将偏微分方程离散化的基础。

可以使用等距网格或非等距网格，具体取决于问题的特点。

2.离散化：根据偏微分方程的类型和边界条件，将偏微分方程的导数转换为离散的差分或有限差分格式。

常用的数值离散化方法有前向差分，后向差分和中心差分等。

3.初值条件：根据问题的初始状态，确定在初始时间步骤上网格点的值。

常用的方法是根据问题的初始条件进行数值插值。

4.边界条件：确定在边界网格点上的值。

根据问题的边界条件，可以采用数值插值法或手动设置边界值。

5. 迭代求解：根据离散化的差分方程，通过迭代方法求解离散化的方程组。

常用的迭代方法有Jacobi方法，Gauss-Seidel方法，SOR方法等。

6.收敛性判断：根据设定的收敛准则，判断数值解是否达到了预期的精度。

通常可以通过比较相邻两次迭代的差异来判断收敛性。

7.后处理：根据求解得到的数值解，计算并绘制出感兴趣的物理量。

还可以评估数值方法的误差和稳定性，并进行必要的修正。

8.参数选择：在数值解的迭代过程中，可能需要选择合适的参数，如网格大小和时间步长等。

这需要根据问题的特性和数值方法的准则进行选择。

9.优化和改进：根据数值解的结果和收敛性，可以对数值方法进行改进和优化。

可能需要调整离散化方法，调整网格布局或改进迭代算法。

总之，偏微分方程的数值解流程是一个迭代过程，通过将偏微分方程离散化为差分方程，并进行迭代求解和收敛性判断，获得问题的数值解。

这个过程需要认真的数值计算和对问题的物理背景知识的深刻理解。

偏微分方程数值解法的MATLAB源码

[原创]偏微分方程数值解法的MATLAB源码【更新完毕】说明:由于偏微分的程序都比较长,比其她的算法稍复杂一些,所以另开一贴,专门上传偏微分的程序谢谢大家的支持！其她的数值算法见:、、//Announce/Announce、asp?BoardID=209&id=82450041、古典显式格式求解抛物型偏微分方程(一维热传导方程)function [U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C)%古典显式格式求解抛物型偏微分方程%[U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C)%%方程:u_t=C*u_xx 0 <= x <= uX,0 <= t <= uT%初值条件:u(x,0)=phi(x)%边值条件:u(0,t)=psi1(t), u(uX,t)=psi2(t)%%输出参数:U -解矩阵,第一行表示初值,第一列与最后一列表示边值,第二行表示第2层……% x -空间变量% t -时间变量%输入参数:uX -空间变量x的取值上限% uT -时间变量t的取值上限% phi -初值条件,定义为内联函数% psi1 -边值条件,定义为内联函数% psi2 -边值条件,定义为内联函数% M -沿x轴的等分区间数% N -沿t轴的等分区间数% C -系数,默认情况下C=1%%应用举例:%uX=1;uT=0、2;M=15;N=100;C=1;%phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0');%[U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C);%设置参数C的默认值if nargin==7C=1;end%计算步长dx=uX/M;%x的步长dt=uT/N;%t的步长x=(0:M)*dx;t=(0:N)*dt;r=C*dt/dx/dx;%步长比r1=1-2*r;if r > 0、5disp('r > 0、5,不稳定')end%计算初值与边值U=zeros(M+1,N+1);for i=1:M+1U(i,1)=phi(x(i));endfor j=1:N+1U(1,j)=psi1(t(j));U(M+1,j)=psi2(t(j));end%逐层求解for j=1:Nfor i=2:MU(i,j+1)=r*U(i-1,j)+r1*U(i,j)+r*U(i+1,j);endendU=U';%作出图形mesh(x,t,U);title('古典显式格式,一维热传导方程的解的图像') xlabel('空间变量x')ylabel('时间变量t')zlabel('一维热传导方程的解U')return;古典显式格式不稳定情况古典显式格式稳定情况2、古典隐式格式求解抛物型偏微分方程(一维热传导方程)function [U x t]=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C)%古典隐式格式求解抛物型偏微分方程%[U x t]=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C)%%方程:u_t=C*u_xx 0 <= x <= uX,0 <= t <= uT%初值条件:u(x,0)=phi(x)%边值条件:u(0,t)=psi1(t), u(uX,t)=psi2(t)%%输出参数:U -解矩阵,第一行表示初值,第一列与最后一列表示边值,第二行表示第2层……% x -空间变量% t -时间变量%输入参数:uX -空间变量x的取值上限% uT -时间变量t的取值上限% phi -初值条件,定义为内联函数% psi1 -边值条件,定义为内联函数% psi2 -边值条件,定义为内联函数% M -沿x轴的等分区间数% N -沿t轴的等分区间数% C -系数,默认情况下C=1%%应用举例:%uX=1;uT=0、2;M=50;N=50;C=1;%phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0');%[U x t]=PDEParabolicClassicalImplicit(uX,uT,phi,psi1,psi2,M,N,C);%设置参数C的默认值if nargin==7C=1;end%计算步长dx=uX/M;%x的步长dt=uT/N;%t的步长x=(0:M)*dx;t=(0:N)*dt;r=C*dt/dx/dx;%步长比Diag=zeros(1,M-1);%矩阵的对角线元素Low=zeros(1,M-2);%矩阵的下对角线元素Up=zeros(1,M-2);%矩阵的上对角线元素for i=1:M-2Diag(i)=1+2*r;Low(i)=-r;Up(i)=-r;endDiag(M-1)=1+2*r;%计算初值与边值U=zeros(M+1,N+1);for i=1:M+1U(i,1)=phi(x(i));endfor j=1:N+1U(1,j)=psi1(t(j));U(M+1,j)=psi2(t(j));end%逐层求解,需要使用追赶法(调用函数EqtsForwardAndBackward)for j=1:Nb1=zeros(M-1,1);b1(1)=r*U(1,j+1);b1(M-1)=r*U(M+1,j+1);b=U(2:M,j)+b1;U(2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b);endU=U';%作出图形mesh(x,t,U);title('古典隐式格式,一维热传导方程的解的图像')xlabel('空间变量x')ylabel('时间变量t')zlabel('一维热传导方程的解U')return;此算法需要使用追赶法求解三对角线性方程组,这个算法在上一篇帖子中已经给出,为了方便,再给出来追赶法解三对角线性方程组function x=EqtsForwardAndBackward(L,D,U,b)%追赶法求解三对角线性方程组Ax=b%x=EqtsForwardAndBackward(L,D,U,b)%x:三对角线性方程组的解%L:三对角矩阵的下对角线,行向量%D:三对角矩阵的对角线,行向量%U:三对角矩阵的上对角线,行向量%b:线性方程组Ax=b中的b,列向量%%应用举例:%L=[-1 -2 -3];D=[2 3 4 5];U=[-1 -2 -3];b=[6 1 -2 1]';%x=EqtsForwardAndBackward(L,D,U,b)%检查参数的输入就是否正确n=length(D);m=length(b);n1=length(L);n2=length(U);if n-n1 ~= 1 || n-n2 ~= 1 || n ~= m disp('输入参数有误！')x=' ';return;end%追的过程for i=2:nL(i-1)=L(i-1)/D(i-1);D(i)=D(i)-L(i-1)*U(i-1);endx=zeros(n,1);x(1)=b(1);for i=2:nx(i)=b(i)-L(i-1)*x(i-1);end%赶的过程x(n)=x(n)/D(n);for i=n-1:-1:1x(i)=(x(i)-U(i)*x(i+1))/D(i);endreturn;古典隐式格式在以后的程序中,我们都取C=1,不再作为一个输入参数处理3、Crank-Nicolson隐式格式求解抛物型偏微分方程需要调用追赶法的程序function [U x t]=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N)%Crank-Nicolson隐式格式求解抛物型偏微分方程%[U x t]=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N)%%方程:u_t=u_xx 0 <= x <= uX,0 <= t <= uT%初值条件:u(x,0)=phi(x)%边值条件:u(0,t)=psi1(t), u(uX,t)=psi2(t)%%输出参数:U -解矩阵,第一行表示初值,第一列与最后一列表示边值,第二行表示第2层……% x -空间变量% t -时间变量%输入参数:uX -空间变量x的取值上限% uT -时间变量t的取值上限% phi -初值条件,定义为内联函数% psi1 -边值条件,定义为内联函数% psi2 -边值条件,定义为内联函数% M -沿x轴的等分区间数% N -沿t轴的等分区间数%%应用举例:%uX=1;uT=0、2;M=50;N=50;%phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0');%[U x t]=PDEParabolicCN(uX,uT,phi,psi1,psi2,M,N);%计算步长dx=uX/M;%x的步长dt=uT/N;%t的步长x=(0:M)*dx;t=(0:N)*dt;r=dt/dx/dx;%步长比Diag=zeros(1,M-1);%矩阵的对角线元素Low=zeros(1,M-2);%矩阵的下对角线元素Up=zeros(1,M-2);%矩阵的上对角线元素for i=1:M-2Diag(i)=1+r;Low(i)=-r/2;Up(i)=-r/2;endDiag(M-1)=1+r;%计算初值与边值U=zeros(M+1,N+1);for i=1:M+1U(i,1)=phi(x(i));endfor j=1:N+1U(1,j)=psi1(t(j));U(M+1,j)=psi2(t(j));endB=zeros(M-1,M-1);for i=1:M-2B(i,i)=1-r;B(i,i+1)=r/2;B(i+1,i)=r/2;endB(M-1,M-1)=1-r;%逐层求解,需要使用追赶法(调用函数EqtsForwardAndBackward) for j=1:Nb1=zeros(M-1,1);b1(1)=r*(U(1,j+1)+U(1,j))/2;b1(M-1)=r*(U(M+1,j+1)+U(M+1,j))/2;b=B*U(2:M,j)+b1;U(2:M,j+1)=EqtsForwardAndBackward(Low,Diag,Up,b);endU=U';%作出图形mesh(x,t,U);title('Crank-Nicolson隐式格式,一维热传导方程的解的图像')xlabel('空间变量x')ylabel('时间变量t')zlabel('一维热传导方程的解U')return;Crank-Nicolson隐式格式4、正方形区域Laplace方程Diriclet问题的求解需要调用Jacobi迭代法与Guass-Seidel迭代法求解线性方程组function [U x y]=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,type) %正方形区域Laplace方程的Diriclet边值问题的差分求解%此程序需要调用Jacobi迭代法或者Guass-Seidel迭代法求解线性方程组%[U x y]=PDEEllipseSquareLaplaceDirichlet(ub,phi1,phi2,psi1,psi2,M,type)%%方程:u_xx+u_yy=0 0<=x,y<=ub%边值条件:u(0,y)=phi1(y)% u(ub,y)=phi2(y)% u(x,0)=psi1(x)% u(x,ub)=psi2(x)%%输出参数:U -解矩阵,第一行表示y=0时的值,第二行表示第y=h时的值……% x -横坐标% y -纵坐标%输入参数:ub -变量边界值的上限% phi1,phi2,psi1,psi2 -边界函数,定义为内联函数% M -横纵坐标的等分区间数% type -求解差分方程的迭代格式,若type='Jacobi',采用Jacobi迭代格式% 若type='GS',采用Guass-Seidel迭代格式。

偏微分方程式(Partial

Y.M. Hu, Assistant Professor, Department of Applied Physics, National University of Kaohsiung
Separation of Variables : Use of Fourier Series
2u T 2u c2 2u
Y.M. Hu, Assistant Professor, Department of Applied Physics, National University of Kaohsiung
Elliptic (Diffusion, Equilibrium Problems)
Y.M. Hu, Assistant Professor, Department of Applied Physics, National University of Kaohsiung
G(t) 0 F(0) F(L) 0
F'' (x) kF(x) 0
For k = 0
F ax b 0x 0 0
X
For positive k = μ2 F Aex Be x 0ex 0ex 0 X
For negative k = -p2 F A cospx Bsin px
Hyperbolic (propagation)
Y.M. Hu, Assistant Professor, Department of Applied Physics, National University of Kaohsiung
Parabolic (Time- or space- marching) 時間或空間步推
derivatives:

python求解偏微分方程

python求解偏微分方程偏微分方程（Partial Differential Equations, PDE）是研究连续介质中的许多物理现象所必需的重要数学工具。

PDE 涉及了空间、时间和其它同步变量之间的关系，因此对于有限元分析（FEM）和流体力学等领域来说，具有极为重要的应用价值。

下面我们将简单介绍使用 Python 求解偏微分方程的基本方法。

1. 引入库在 Python 中，我们可以使用 SciPy 和 NumPy 库来处理偏微分方程。

其中，NumPy 用于数值计算，而 SciPy 则提供了一些特定的算法，包括线性方程组求解、优化、数值积分和微分方程等。

因此，我们需要在程序中引入这两个库：```pythonimport numpy as npfrom scipy import sparsefrom scipy.sparse.linalg import spsolve```2. 构建矩阵在求解偏微分方程时，我们通常需要构建雅可比矩阵。

这里举一个简单的例子，设有一个一维热传导方程：$$ \frac{\partial^2 u}{\partial x^2} = f(x) $$其中，$u$ 是未知函数，$f(x)$ 是给定函数。

为了求解这个方程，我们可以按照离散化的方法来处理。

我们将区间 $[0,1]$ 分成 $n$ 个小区间，即 $x_0 = 0$，$x_n= 1$，$x_i = ih$，$h = 1/n$。

因此，$u(x_i)$ 可以用 $u_i$ 来表示。

我们将前式中的二阶导数离散化，得到如下近似式：$$ \frac{u_{i+1} - 2u_i + u_{i-1}}{h^2} = f_i $$解出 $u_i$，得到：$$u_i = \frac{1}{h^2} \left(u_{i+1} + u_{i-1} - h^2f_i\right)$$这样，我们就可以得到一个线性方程组：$$A\mathbf{u} = \mathbf{f}$$其中，$\mathbf{u}$ 是 $[u_1, u_2, ..., u_n]$ 的列向量，$\mathbf{f}$ 是 $[f_1, f_2, ..., f_n]$ 的列向量。

ADI(交替隐式求求解抛物型偏微分方程)程序

% 精确解: u(t,x,y) = exp( 1/2*(x+y) - t ) %
% 右端函数: f(t,x,y) = 1/2*exp( 1/2*(x+y) - t ) %
% 初始条件: g1 = exp( 1/2*(x+y) ) %
h = (b-a)/n; dt=h;
h1 = h*h;
r=dt/h1;
x=a:h:b; y=c:h:d;
%-- Initial condition:
t = 0;
for i=1:m+1,
for j=1:m+1,
u1(i,j) = uexact(t,x(i),y(j));
% %
% Files needed for the test: %
% %
% %
% Example of ADI Method for 2D heat equation %
% %
% u_tt = u_{xx} + u_{yy} + f(x,t) %
% a<x<b;c<y<d %
% Test problme: %
end %-- Finished with the loop in time
%----------- Data analysis ----------------------------------
for i=1:m+1,
for j=1:n+1,
ue(i,j) = uexact(tfinal,x(i),y(j));
e = max(max(abs(u1-ue))) % The infinity error.

偏微分方程的数值解方法及源程序

-240-第二十章偏微分方程的数值解自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。

这些规律的定量表述一般地呈现为关于含有未知函数及其导数的方程。

我们将只含有未知多元函数及其偏导数的方程，称之为偏微分方程。

方程中出现的未知函数偏导数的最高阶数称为偏微分方程的阶。

如果方程中对于未知函数和它的所有偏导数都是线性的，这样的方程称为线性偏微分方程，否则称它为非线性偏微分方程。

初始条件和边界条件称为定解条件，未附加定解条件的偏微分方程称为泛定方程。

对于一个具体的问题，定解条件与泛定方程总是同时提出。

定解条件与泛定方程作为一个整体，称为定解问题。

§1 偏微分方程的定解问题各种物理性质的定常（即不随时间变化）过程，都可用椭圆型方程来描述。

其最典型、最简单的形式是泊松(Poisson)方程),(2222y x f y ux u u =∂∂+∂∂=Δ （1）特别地，当0),(≡y x f 时，即为拉普拉斯(Laplace)方程，又称为调和方程02222=∂∂+∂∂=Δyux u u （2）带有稳定热源或内部无热源的稳定温度场的温度分布，不可压缩流体的稳定无旋流动及静电场的电势等均满足这类方程。

Poisson 方程的第一边值问题为⎪⎩⎪⎨⎧Ω∂=Γ=Ω∈=∂∂+∂∂Γ∈),(|),(),(),(),(2222y x y x u y x y x f y uxu y x ϕ （3）其中Ω为以Γ为边界的有界区域，Γ为分段光滑曲线，ΓΩU 称为定解区域，),(),,(y x y x f ϕ分别为ΓΩ,上的已知连续函数。

第二类和第三类边界条件可统一表示成),(),(y x u n u y x ϕα=⎟⎠⎞⎜⎝⎛+∂∂Γ∈ （4）其中n 为边界Γ的外法线方向。

当0=α时为第二类边界条件，0≠α时为第三类边界条件。

在研究热传导过程，气体扩散现象及电磁场的传播等随时间变化的非定常物理问题时，常常会遇到抛物型方程。

偏微分方程数值解程序作业

偏微分方程数值解程序作业二维热传导#include "iostream"#include "math.h"#include "stdio.h"#define TT 0.1//要求的t的值#define PI 3.1415926535using namespace std;#define N 11//跳点全局变量double dot_sheet_b[N][N];double dot_sheet_t[N][N];//显示全局变量double dot_sheet_old[N][N];double dot_sheet_new[N][N];double lamda=0.5;//lamda=a*t/powf(h,2) double h1=0.1,h2=0.1;void get_next_dot(int i,int j);void initial(void);void initial(void);void initial1(void);void get_next_2_layer(void);//每次得到两层void get_next_1_layer(void);//显示一次获得一层void main(void){int i,j;//行列计数使用initial();get_next_2_layer();cout<<"output_line"<<endl;for( i=1;i<N-1;i++)//输出跳点格式的（半隐式）{for ( j=1;j<N-1;j++)printf("%.3f|",dot_sheet_b[i][j]);cout<<endl;}printf("**************隐式结果***************\n");initial1();get_next_1_layer();for( i=1;i<N-1;i++)//输出显示方式的{for ( j=1;j<N-1;j++)printf("%.3f|",dot_sheet_new[i][j]);cout<<endl;}printf("***************显示结果*******************\n"); }void get_next_1_layer(void)//显示一次获得一层{double delt_x2,delt_y2;int i,j;for( i=1;i<N-1;i++){for ( j=1;j<N-1;j++){delt_x2=(dot_sheet_old[i+1][j]-2*dot_sheet_old[i][j]+dot_sheet_old[i-1][j]); delt_y2=(dot_sheet_old[i][j+1]-2*dot_sheet_old[i][j]+dot_sheet_old[i][j-1]); dot_sheet_new[i][j]=dot_sheet_old[i][j]+lamda*(delt_x2+delt_y2);}}}void get_next_2_layer(void)//跳点格式一次获得两层{double delt_x2,delt_y2;int i,j;for ( i=1;i<N-1;i++)////奇数，显示方式{for(j=(i+1)%2+1;j<N-1;j+=2){delt_x2=(dot_sheet_t[i+1][j]-2*dot_sheet_t[i][j]+dot_sheet_t[i-1][j]);delt_y2=(dot_sheet_t[i][j+1]-2*dot_sheet_t[i][j]+dot_sheet_t[i][j-1]);dot_sheet_b[i][j]=dot_sheet_t[i][j]+lamda*(delt_x2+delt_y2);// printf("%.3f|",dot_sheet_b[i][j]);}// cout<<endl;}for ( i=1;i<N-1;i++)////偶数，显示方式{for(j=(i)%2+1;j<N-1;j+=2){delt_x2=(dot_sheet_t[i+1][j]-2*dot_sheet_t[i][j]+dot_sheet_t[i-1][j]);delt_y2=(dot_sheet_t[i][j+1]-2*dot_sheet_t[i][j]+dot_sheet_t[i][j-1]);dot_sheet_b[i][j]=lamda*(delt_x2+delt_y2)+dot_sheet_t[i][j];}// cout<<endl;}cout<<"***********************************************************"<<endl;for (i=1;i<N-1;i++)//第2层{for(j=1;j<N-1;j++){dot_sheet_t[i][j]=2*dot_sheet_b[i][j]-dot_sheet_t[i][j];// printf("%.4f|",dot_sheet_t[i][j]);//第二层输出屏蔽，暂时只输出第一层}// cout<<endl;}}void initial(void)//跳点格式初始化{for(int i=0;i<N;i++)dot_sheet_t[i][0]=dot_sheet_t[i][N-1]=0;for (int j=0;j<N;j++)dot_sheet_t[0][j]=dot_sheet_t[N-1][j]=0;for (i=1;i<N-1;i++){for ( j=1;j<N-1;j++){dot_sheet_b[i][j]=dot_sheet_t[i][j]=sin((h1*i+h2*j)*PI);}cout<<endl;}printf("************初始化底层**********************\n"); }void initial1(void)//显示初始化{for(int i=0;i<N;i++)dot_sheet_old[i][0]=dot_sheet_old[i][N-1]=0;for (int j=0;j<N;j++)dot_sheet_old[0][j]=dot_sheet_old[N-1][j]=0;for (i=1;i<N-1;i++){for ( j=1;j<N-1;j++){dot_sheet_old[i][j]=sin((h1*i+h2*j)*PI);}}printf("**************初始化底层1*******************\n"); }一维热传导-----向前差分#include "iostream"#include "math.h"#include "stdio.h"#define TT 0.1//要求的t的值#define PI 3.1415926535using namespace std;void main(void){double r=0.1;//0.5double step_h=0.1;double step_t;//初始化无任何意义int row1,cln1,row_max;step_t=r*(step_h*step_h);// printf("row_max=%f\n",step_t);row_max=TT/step_t+1;//认为是第rowmax行printf("row_max=%d\n",row_max);double matrix_u[20][11];//一行共11个点//以下为带入边界条件for (cln1=0;cln1<=10;cln1++)matrix_u[0][cln1]=sin(PI*step_h*cln1);for(row1=0;row1<=row_max;row1++)matrix_u[row1][0]=matrix_u[row1][10]=0;for (row1=1;row1<=row_max;row1++)for(cln1=1;cln1<=9;cln1++)matrix_u[row1][cln1]=r*(matrix_u[row1+1][cln1-1]+matrix_u[row1-1][cln1-1])+(1-r)*(matrix_ u[row1][cln1-1]);for(cln1=0;cln1<=10;cln1++){printf("matrix_u[%d][%d] =%.5f\n",row1-1,cln1,matrix_u[row1-1][cln1]);printf("the exact value=%.5f\n",exp(-PI*PI*(row_max-1)*step_t)*sin(PI*step_h*cln1));}}向后差分格式#include "iostream"#include "math.h"#include "stdio.h"#include "malloc.h"#define TT 0.1//要求的t的值#define PI 3.1415926535using namespace std;void zhuigan(double *A,double *B,double *C,double *D,int n)//n表示矩阵的维数，就是行数。

偏微分方程的数值解方法及源程序

如果偏微分方程定解问题的解存在，唯一且连续依赖于定解数据（即出现在方程和定解条件中的已知函数），则此定解问题是适定的。可以证明，上面所举各种定解问题都是适定的。 §2 偏微分方程的差分解法差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。如果差分格式有解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为原问题的近似解（数值解）。因此，用差分方法求偏微分方程定解问题一般需要解决以下问题：（i）选取网格；（ii）对微分方程及定解条件选择差分近似，列出差分格式；（iii）求解差分格式；（iv）讨论差分格式解对于微分方程解的收敛性及误差估计。下面我们只对偏微分方程的差分解法作一简要的介绍。 2.1 椭圆型方程第一边值问题的差分解法以 Poisson 方程（1）为基本模型讨论第一边值问题的差分方法。考虑 Poisson 方程的第一边值问题（3）
（8）
∂u ∂u +a =0 ∂t ∂x
物理中常见的一维振动与波动问题可用二阶波动方程
2 ∂ 2u 2 ∂ u =a ∂t 2 ∂x 2
（9）
（10）
描述，它是双曲型方程的典型形式。方程（10）的初值问题为
2 ⎧ ∂ 2u 2 ∂ u ⎪ 2 =a ∂t ∂x 2 ⎪ ⎪ ⎨u ( x,0) = ϕ ( x) ⎪ ∂u ⎪ = φ ( x) ⎪ ⎩ ∂t t =0

偏微分方程程序附件

程序一、计算下面双曲方程的近似解解法：1、一阶迎风格式程序：function upwind1(h,lamda) %迎风格式一阶精度显示格式 x=h:h:10-h;tao=lamda*h;t=tao:tao:5-tao;[X,T]=meshgrid(x,t);i=1:10/h;j=1:5/tao;u(i,j)=0;for n=2:5/taofor m=2:10/hu(m,n)=u(m,n-1)-lamda*(u(m,n-1)-u(m-1,n-1))+2*tao*(n-1)*tao*sin(m*h)+(n-1)^2*tao^3*cos(m *h); %迎风格式endendfor n=1:5/tao-1for m=1:10/h-1U(m,n)=u(m+1,n+1);endendU=U';figure(1);grid on;mesh(X,T,U);xlabel('x 轴');ylabel('t 轴');zlabel('u 轴');title('迎风格式差分曲面图');%计算误差⎪⎪⎩⎪⎪⎨⎧==<<<<+=∂∂+∂∂== 0| 0|105,00,cos sin 2 002x t u u x t x t x t x u t ujj=1:5/tao-1;uu(ii,jj)=0;for nn=1:5/tao-1for mm=1:10/h-1uu(mm,nn)=(nn*tao)^2*sin(mm*h);endendUU=uu';error= U-UU;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('迎风格式误差曲面图');figure(3);mesh(X,T,UU);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('古典解');2、Lax-Wendroff格式function lax_w1(h,lamda) %Lax-Wendroff 二阶精度显示格式x=h:h:10-h;tao=lamda*h;t=tao:tao:5-tao;[X,T]=meshgrid(x,t);i=1:10/h;j=1:5/tao;u(i,j)=0;for n=2:5/taofor m=2:10/hif m==10/hu(m,n)=u(m,n-1)-lamda*(u(m,n-1)-u(m-1,n-1))+2*tao*(n-1)*tao*sin(m*h)...+(n-1)^2*tao^3*cos(m*h); %用迎风格式处理边界问题elseu(m,n)=u(m,n-1)-0.5*lamda*(u(m+1,n-1)-u(m-1,n-1))+2*tao*(n-1)*tao*sin(m*h)...+tao*((n-1)*tao)^2*cos(m*h)+0.5*lamda^2*(u(m+1,n-1)-2*u(m,n-1)+u(m-1,n-1))... +tao^2*(n-1)*tao*(sin(m*h)-cos(m*h))+0.5*tao^2*(n-1)^2*tao^2*(sin(m*h)+cos(m*h)); %Lax-Wendroff差分格式endendfor n=1:5/tao-1for m=1:10/h-1U(m,n)=u(m+1,n+1);endendU=U';figure(1);grid on;mesh(X,T,U);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('Lax-Wendroff格式差分曲面图');%计算误差ii=1:10/h-1;jj=1:5/tao-1;uu(ii,jj)=0;for nn=1:5/tao-1for mm=1:10/h-1uu(mm,nn)=(nn*tao)^2*sin(mm*h);endendUU=uu';error=U-UU;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('Lax-Wendroff格式误差曲面图');3、隐式中心格式function implicit_1(h,lamda) %隐式格式，精度为：时间一阶，空间二阶；无条件稳定x=h:h:10-h;tao=lamda*h;t=tao:tao:5-tao;[X,T]=meshgrid(x,t);time1=5/tao; %时间为1时的时间层网格点space1=10/h; %空间为1时的时间层网格点j=1:time1;u(i,j)=0;r=lamda/2;%使用追赶法求每一时间层的节点值k=1:space1-1; %初始化系数矩阵a(k,k)=0;a(1,1)=1;a(1,2)=r;for k=2:space1-2 %定义系数矩阵a(k,k-1)=-r;a(k,k)=1;a(k,k+1)=r;enda(space1-1,space1-1)=1;p=1:space1-1;q=1:time1-1;Z(p,q)=0; %矩阵a*x(i)=Z(i);for n=2:time1for m=2:space1-1Z(m-1,n-1)=u(m,n-1)+tao*(2*(n-1)*tao*sin((m-1)*h)+...n^2*tao^2*cos(m*h));endZ(space1-1,n-1)=u(space1,n-1)-lamda*(u(space1,n-1)-u(space1-1,n-1))...+2*tao*(n-1)*tao*sin(space1*h)+...(n-1)^2*tao^3*cos(space1*h); %边界离散，时间一阶，空间一阶E=chase(a,Z(:,n-1)); %调用追赶法程序，程序在下面for v=1:space1-1u(v+1,n)=E(v);endendfor n=1:5/tao-1for m=1:10/h-1U(m,n)=u(m+1,n+1);endendU=U';figure(1);grid on;mesh(X,T,U);ylabel('t轴');zlabel('u轴');title('隐式格式差分曲面图');%计算误差ii=1:10/h-1;jj=1:5/tao-1;uu(ii,jj)=0;for nn=1:5/tao-1for mm=1:10/h-1uu(mm,nn)=(nn*tao)^2*sin(mm*h);endendUU=uu';error=UU-U;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('隐式格式误差曲面图');调用追赶法接三对角矩阵函数chase.m为：function X=chase(A,d)%追赶法只用来解三对角方程组，其中A为方程的系数，d为右端项%a为-1对角线上元素%b为主对角线上的元素%c为+1对角线上的元素a=diag(A,-1);b=diag(A);c=diag(A,1);n=length(b);U(1)=b(1);y(1)=d(1);for i=2:nL(i-1)=a(i-1)/U(i-1);U(i)=b(i)-c(i-1)*L(i-1);y(i)=d(i)-L(i-1)*y(i-1);endX(n)=y(n)/U(n);for i=n-1:-1:1X(i)=(y(i)-c(i)*X(i+1))/U(i);end二、计算下面双曲方程的近似解解法：1、二阶显示格式程序：function explicit_2(h,lamda) %波动方程显示格式，二阶精度，lamda<1时稳定 tao=h*lamda;t=tao:tao:1-tao;x=h:h:1-h;time1=1/tao+1; %时间为1时的时间层网格点space1=1/h+1; %空间为1时的时间层网格点[X,T]=meshgrid(x,t);i=1:space1;j=1:time1;u(i,j)=0;for n=2:time1 %边界条件离散u(space1,n)=sin((n-1)*tao);endfor m=2:1/h %初始条件离散二阶精度u(m,2)=u(m,1)+tao*h*(m-1)+0.5*lamda^2*(u(m+1,1)-2*u(m,1)+u(m-1,1));endfor n=3:1/taofor m=2:space1-1u(m,n)=2*u(m,n-1)-u(m,n-2)+lamda^2*(u(m+1,n-1)-2*u(m,n-1)+u(m-1,n-1))+... tao^2*(((n-1)*tao)^2-((m-1)*h)^2)*sin((n-1)*tao*(m-1)*h); %波动方程显示格式 end end⎪⎪⎪⎩⎪⎪⎪⎨⎧≤<==<<=∂∂=<<<<-=∂∂∂∂ 10 ,sin t)u(1,, 0),0(10 ,,0)(t , 0)0,(11,00 ),sin()(- 222222t t t u x x x u x u x t xt x t x u t ufor n=1:time1-2for m=1:space1-2U(m,n)=u(m+1,n+1);endendU=U';figure(1);grid on;mesh(X,T,U);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('波动方程显示格式差分曲面图');%计算古典解ii=1:space1-2;jj=1:time1-2;uu(ii,jj)=0;for nn=1:time1-2for mm=1:space1-2uu(mm,nn)=sin(mm*h*nn*tao);endendUU=uu';error= U-UU;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('波动方程显示格式误差曲面图');figure(3);mesh(X,T,UU);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('波动方程古典解');2、隐式二阶格式程序：function implicit_2(h,lamda) %波动方程隐式格式，二阶精度，无条件稳定tao=h*lamda;t=tao:tao:1-tao;x=h:h:1-h;time1=1/tao+1; %时间为1时的时间层网格点space1=1/h+1; %空间为1时的时间层网格点[X,T]=meshgrid(x,t);i=1:space1;j=1:time1;u(i,j)=0; %离散化的网格点v(i,j)=0; %定义关于时间的一阶导数for i=1:space1;v(i,1)=(i-1)*h; %时间的一阶导数初始值离散化endfor i=1:space1-1for j=1:time1w(i,j)=0; %定义关于空间的一阶导数endendfor i=1:space1-1w(i,1)=0; %空间的一阶导数初始值离散化endfor n=2:time1 %边界条件离散u(space1,n)=sin((n-1)*tao);end%使用追赶法求每一时间层v的节点值,需要定义系数矩阵a和非其次项Z r=-0.25*lamda^2;m=1+0.5*lamda^2;k=1:space1-2; %初始化系数矩阵a(k,k)=0;a(1,1)=m;a(1,2)=r;for k=2:space1-3 %定义系数矩阵a(k,k-1)=r;a(k,k)=m;a(k,k+1)=r;enda(space1-2,space1-3)=r;a(space1-2,space1-2)=m;%定义非齐次项p=1:space1-2;q=1:time1-1;Z(p,q)=0; %矩阵a*x(i)=Z(i);%使用隐式格式计算for jj=2:time1 %从第二时间层开始计算v(1,jj)=0; %v左边界值离散v(space1,jj)=(u(space1,jj)-u(space1,jj-1))/tao; %v右边界值离散%计算非齐次项for m=1:space1-3Z(m,jj-1)=v(m+1,jj-1)+lamda*(w(m+1,jj-1)-w(m,jj-1))-...r*(v(m+2,jj-1)-2*v(m+1,jj-1)+v(m,jj-1));endZ(space1-2,jj-1)=-r*v(space1,jj)+v(space1-1,jj-1)+lamda*...(w(space1-1,jj-1)-w(space1-2,jj-1))-...r*(v(space1,jj-1)-2*v(space1-1,jj-1)+v(space1-2,jj-1));E=chase(a,Z(:,jj-1));for kk=1:space1-2v(kk+1,jj)=E(kk);end%计算w(:,jj)和u(:,jj)for ii=1:space1-1w(ii,jj)=0.5*lamda*(v(ii+1,jj)-v(ii,jj))+w(ii,jj-1)+...0.5*lamda*(v(ii+1,jj-1)-v(ii,jj-1));endfor g=2:space1-1u(g,jj)=u(g,jj-1)+tao*v(g,jj);endendfor n=1:time1-2for m=1:space1-2U(m,n)=u(m+1,n+1);endendU=U';figure(1);grid on;mesh(X,T,U);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('波动方程显示格式差分曲面图');%计算古典解i1=1:space1-2;jj=1:time1-2;uu(i1,jj)=0;for nn=1:time1-2for mm=1:space1-2uu(mm,nn)=sin(mm*h*nn*tao);endendUU=uu';error= U-UU;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('波动方程显示格式误差曲面图');三、计算下面抛物方程的近似解1、向前差分格式程序：function diffusion_1(h,lamda,time)%扩散方程向前差分格式，时间一阶精度，空间二阶精度，lamda<=0,5稳定 %h 为空间步长，lamda 为网格比，t 为演化时间%h=0.1;time=1;lamda=0.5;tao=h^2*lamda;t=0:tao:time;x=0:h:1;[X,T]=meshgrid(x,t);time1=time/tao+2; %时间为1时的时间层网格点space1=1/h+1; %空间为1时的时间层网格点i=1:space1;j=1:time1;u(i,j)=0; %离散化的网格点wucha(1,i)=0;for i=1:space1-1 %边界条件初始化u(i,1)=sin(pi*(i-1)*h);endfor n=2:time1for m=2:space1-1u(m,n)=(1-2*lamda)*u(m,n-1)+lamda*u(m-1,n-1)+lamda*u(m+1,n-1); endendu=u';figure(1);grid on;mesh(X,T,u);xlabel('x 轴');ylabel('t 轴');zlabel('u 轴');title('向前差分格式曲面图');%计算误差ii=1:space1; jj=1:time1;⎪⎪⎪⎩⎪⎪⎪⎨⎧<==<<=<<<∂∂=∂∂ 0 ,0t)u(1,),0(10 sin )0,(1,00 22t t u x x x u x t x u t u πuu=exp(-pi^2.*T).*sin(pi.*X);error=uu-u;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('向前差分格式误差曲面图');%plot(error);figure(3);mesh(X,T,uu);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('古典解');wucha=error(22,:)figure(4)plot(wucha,'ro-');xlabel('t=0.1s，x轴');ylabel('误差轴');title('0.1秒时差分格式解与精确解的误差');2、向后差分格式程序function diffusion_2(h,lamda,time)%扩散方程向后差分格式，时间一阶精度，空间二阶精度，无条件稳定%h为空间步长，lamda为网格比，t为演化时间%h=0.1;time=1;lamda=0.5;tao=h^2*lamda;t=0:tao:time;x=0:h:1;[X,T]=meshgrid(x,t);time1=time/tao+2; %时间为1时的时间层网格点space1=1/h+1; %空间为1时的时间层网格点i=1:space1;j=1:time1;u(i,j)=0; %离散化的网格点wucha(1,i)=0;for i=1:space1-1 %边界条件初始化u(i,1)=sin(pi*(i-1)*h);end%使用追赶法求每一时间层u的节点值,需要定义系数矩阵a和非其次项Zm=1+2*lamda;kk=1:space1-2; %初始化系数矩阵a(kk,kk)=0;a(1,1)=m;a(1,2)=r;for k=2:space1-3 %定义系数矩阵a(k,k-1)=r;a(k,k)=m;a(k,k+1)=r;enda(space1-2,space1-3)=r;a(space1-2,space1-2)=m;%定义非齐次项p=1:space1-2;q=1:time1-1;Z(p,q)=0; %矩阵a*x(i)=Z(i); %使用隐式格式计算for jj=2:time1 %从第二时间层开始计算%计算非齐次项for m=1:space1-2Z(m,jj-1)=u(m+1,jj-1);end%追赶法E=chase(a,Z(:,jj-1));for kk=1:space1-2u(kk+1,jj)=E(kk);endendu=u';figure(1);grid on;mesh(X,T,u);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('向后差分格式曲面图');%计算误差uu=exp(-pi^2.*T).*sin(pi.*X);error=uu-u;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel('x轴');zlabel('u轴');title('向后差分格式误差曲面图');wucha=error(22,:)figure(3)plot(wucha,'ro-');xlabel('t=0.1s，x轴');ylabel('误差轴');title('0.1秒时向后差分格式解与精确解的误差');3、Crank-Nicolson格式程序function diffusion_3(h,lamda,time)%扩散方程Crank-Nicolson差分格式，时间二阶精度，空间二阶精度，无条件稳定%h为空间步长，lamda为网格比，t为演化时间%h=0.1;time=1;lamda=0.5;tao=h^2*lamda;t=0:tao:time;x=0:h:1;[X,T]=meshgrid(x,t);time1=time/tao+2; %时间为1时的时间层网格点space1=1/h+1; %空间为1时的时间层网格点i=1:space1;j=1:time1;u(i,j)=0; %离散化的网格点wucha(1,i)=0;for i=1:space1-1 %边界条件初始化u(i,1)=sin(pi*(i-1)*h);end%使用追赶法求每一时间层u的节点值,需要定义系数矩阵a和非其次项Zr=-0.5*lamda;r1=-r;m=1+lamda;m1=1-lamda;kk=1:space1-2; %初始化系数矩阵a(kk,kk)=0;a(1,1)=m;a(1,2)=r;for k=2:space1-3 %定义系数矩阵a(k,k-1)=r;a(k,k)=m;a(k,k+1)=r;a(space1-2,space1-3)=r;a(space1-2,space1-2)=m;%定义非齐次项p=1:space1-2;q=1:time1-1;Z(p,q)=0; %矩阵a*x(i)=Z(i);%使用隐式格式计算for jj=2:time1 %从第二时间层开始计算%计算非齐次项for m=1:space1-2Z(m,jj-1)=r1*u(m,jj-1)+m1*u(m+1,jj-1)+r1*u(m+2,jj-1);end%追赶法E=chase(a,Z(:,jj-1));for kk=1:space1-2u(kk+1,jj)=E(kk);endendu=u';figure(1);grid on;mesh(X,T,u);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('Crank-Nicolson差分格式曲面图');%计算误差uu=exp(-pi^2.*T).*sin(pi.*X);error=uu-u;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('Crank-Nicolson差分格式误差曲面图');wucha=error(22,:)figure(3)plot(wucha,'ro-');xlabel('t=0.1s，x轴');ylabel('误差轴');title('0.1秒时Crank-Nicolson差分格式解与精确解的误差');4、Du Fort-Frankel格式程序function diffusion_4(h,lamda,time)%扩散方程Du Fort-Frankel差分格式，时间二阶精度，空间二阶精度，无条件稳定%h为空间步长，lamda为网格比，t为演化时间%h=0.1;time=1;lamda=0.5;tao=h^2*lamda;t=0:tao:time;x=0:h:1;[X,T]=meshgrid(x,t);time1=time/tao+2; %时间为1时的时间层网格点space1=1/h+1; %空间为1时的时间层网格点i=1:space1;j=1:time1;u(i,j)=0; %离散化的网格点wucha(1,i)=0;for i=1:space1-1 %边界条件初始化u(i,1)=sin(pi*(i-1)*h);end%使用追赶法求每一时间层u的节点值,需要定义系数矩阵a和非其次项Zr=-0.5*lamda;r1=-r;m=1+lamda;m1=1-lamda;kk=1:space1-2; %初始化系数矩阵a(kk,kk)=0;a(1,1)=m;a(1,2)=r;for k=2:space1-3 %定义系数矩阵a(k,k-1)=r;a(k,k)=m;a(k,k+1)=r;enda(space1-2,space1-3)=r;a(space1-2,space1-2)=m;%定义非齐次项p=1:space1-2;q=1;Z(p,q)=0; %矩阵a*x(i)=Z(i);jj=2; %使用Crank-Nicolson差分格式计算第二时间层%计算非齐次项for m=1:space1-2Z(m,jj-1)=r1*u(m,jj-1)+m1*u(m+1,jj-1)+r1*u(m+2,jj-1);end%追赶法E=chase(a,Z(:,jj-1));for kk=1:space1-2u(kk+1,jj)=E(kk);end%使用Du Fort-Frankel差分格式计算u值d1=(1-2*lamda)/(1+2*lamda);d2=2*lamda/(1+2*lamda);for j1=3:time1for i1=2:space1-1u(i1,j1)=d1*u(i1,j1-2)+d2*u(i1-1,j1-1)+d2*u(i1+1,j1-1);endendu=u';figure(1);grid on;mesh(X,T,u);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('Crank-Nicolson差分格式曲面图');%计算误差uu=exp(-pi^2.*T).*sin(pi.*X);error=uu-u;figure(2);mesh(X,T,error);%mesh(X,T,UU);xlabel('x轴');ylabel('t轴');zlabel('u轴');title('Crank-Nicolson差分格式误差曲面图');wucha=error(22,:)figure(3)plot(wucha,'ro-');xlabel('t=0.1s，x轴');ylabel('误差轴');title('0.1秒时Crank-Nicolson差分格式解与精确解的误差');。

一维抛物线偏微分方程数值解法(3)(附图及matlab程序)培训讲学

一维抛物线偏微分方程数值解法(3)(附图及m a t l a b程序)一维抛物线偏微分方程数值解法（3）上一篇参看一维抛物线偏微分方程数值解法（2）（附图及matlab程序）解一维抛物线型方程(理论书籍可以参看孙志忠：偏微分方程数值解法）Ut-Uxx=0, 0<x<1,0<t<=1（Ut-aUxx=f(x,t),a>0)U(x,0)=e^x, 0<=x<=1,U(0,t)=e^t,U(1,t)=e^(1+t), 0<t<=1精确解为：U(x,t)=e^(x+t);此种方法精度为o(h1^2+h2^2)一：用追赶法解线性方程组（还可以用迭代法解）Matlab程序function [u p e x t]=CN(h1,h2,m,n)%Crank-Nicolson格式差分法解一维抛物线型偏微分方程%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m,n分别为空间，时间网格数%p为精确解，u为数值解,e为误差x=(0:m)*h1+0; x0=(0:m)*h1; %定义x0,t0是为了f(x,t)~=0的情况%t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;syms f;for(i=1:n+1)for(j=1:m+1)f(i,j)=0; %f(i,j)=f(x0(j),t0(i))==0% endendfor(i=1:n+1)u(i,1)=exp(t(i));u(i,m+1)=exp(1+t(i));endfor(i=1:m+1)u(1,i)=exp(x(i));endr=h2/(h1*h1);for(i=1:n) %外循环，先固定每一时间层，每一时间层上解一线性方程组%a(1)=0;b(1)=1+r;c(1)=-r/2;d(1)=r/2*(u(i+1,1)+u(i,1))+h2*f(i,j)...+(1-r)*u(i,2)+r/2*u(i,3);for(k=2:m-2)a(k)=-r/2;b(k)=1+r;c(k)=-r/2;d(k)=h2*f(i,j)+r/2*u(i,k)+(1-r)...*u(i,k+1)+r/2*u(i,k+2);%输入部分系数矩阵，为0的矩阵元素不输入%enda(m-1)=-r/2;b(m-1)=1+r;d(m-1)=h2*f(i,j)+r/2*(u(i,m+1)+u(i+1,m+1)...)+r/2*u(i,m-1)+(1-r)*u(i,m);for(k=1:m-2) %开始解线性方程组消元过程a(k+1)=-a(k+1)/b(k);b(k+1)=b(k+1)+a(k+1)*c(k);d(k+1)=d(k+1)+a(k+1)*d(k);endu(i+1,m)=d(m-1)/b(m-1); %回代过程%for(k=m-2:-1:1)u(i+1,k+1)=(d(k)-c(k)*u(i+1,k+2))/b(k);endendfor(i=1:n+1)for(j=1:m+1)p(i,j)=exp(x(j)+t(i)); %p为精确解e(i,j)=abs(u(i,j)-p(i,j));%e为误差endend[u p e x t]=CN(0.1,0.005,10,200);surf(x,t,e); shading interp;>> xlabel('x');ylabel('t');zlabel('e'); >> title('误差曲面')plot(x,e)plot(t,e)误差较向前欧拉法减小一半但是运行时间较长，约39秒,而前两次运行只需l秒左右；[u p e x t]=CN(0.01,0.01,100,100);运行需三分钟左右,误差比前次提高五倍,运算量也提高五倍[u p e x t]=CN(0.1,0.1,10,10);surf(x,t,e) 运行需要2秒；精度还是挺高的；[u p e x t]=CN(0.1,0.2,10,5);surf(x,t,e)误差还可以接受此种方法精度高，计算量较大二：用迭代法解线性方程组：Matlab程序如下：function [u e p x t k]=CN1(h1,h2,m,n,kmax,ep)% 解抛物线型一维方程 C-N格式（Ut-aUxx=f(x,t),a>0)%用g-s(高斯-赛德尔)迭代法解%kmax为最大迭代次数%m,n为x,t方向的网格数，例如（2-0）/0.01=200;%e为误差，p为精确解syms temp;u=zeros(n+1,m+1);x=0+(0:m)*h1;t=0+(0:n)*h2;for(i=1:n+1)u(i,1)=exp(t(i));u(i,m+1)=exp(1+t(i));endfor(i=1:m+1)u(1,i)=exp(x(i));endfor(i=1:n+1)for(j=1:m+1)f(i,j)=0;endenda=zeros(n,m-1);r=h2/(h1*h1); %此处r=a*h2/(h1*h1）；a=1for(k=1:kmax)for(i=1:n)for(j=2:m)temp=((r/2*u(i,j-1)+(1-r)*u(i,j)+r/2*u(i,... j+1)+h2*f(i,j)+r/2*u(i+1,j-1)+r/2*u(i+1,j+1))/(1+r));a(i+1,j)=(temp-u(i+1,j))*(temp-u(i+1,j));u(i+1,j)=temp;%此处注意是u(i+1,j),,而不是u(i+1,j+1)%endenda(i+1,j)=sqrt(a(i+1,j));if(k>kmax)break;endif(max(max(a))<ep)break;endendfor(i=1:n+1)for(j=1:m+1)p(i,j)=exp(x(j)+t(i));e(i,j)=abs(u(i,j)-p(i,j));endend[u e p x t k]=CN1(0.1,0.005,10,200,10000,1e-10);运行速度：1秒迭代次数k =81surf(x,t,e)第二幅图为三角追赶法解方程作出的图，两者几乎一样；由于迭代法速度很快，所以可以将区间分得更小[u e p x t k]=CN1(0.01,0.01,100,100,10000,1e-12);surf(x,t,e);shading interp; k=6903。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。