偏微分方程解的几道算例(差分、有限元)+含matlab程序
- 格式:pdf
- 大小:443.53 KB
- 文档页数:15
基础知识偏微分方程的定解问题各种物理性质的定常(即不随时间变化)过程,都可用椭圆型方程来描述。
其最典型、最简单的形式是泊松(Poisson)方程),(2222y x f yux u u =∂∂+∂∂=∆ (1)特别地,当 f ( x , y ) ≡ 0 时,即为拉普拉斯(Laplace)方程,又称为调和方程02222=∂∂+∂∂=∆yux u u (2)带有稳定热源或内部无热源的稳定温度场的温度分布,不可压缩流体的稳定无旋流动及静电场的电势等均满足这类方程。
Poisson 方程的第一边值问题为⎪⎩⎪⎨⎧Ω∂=Γ=Ω∈=∂∂+∂∂=∆Γ∈),(),(),(),(),(2222y x y x u y x y x f y ux u u y x ϕ (3) 其 中 Ω 为 以 Γ 为 边 界 的 有 界区 域 , Γ 为 分 段 光 滑 曲 线, Ω U Γ 称 为 定 解区 域 ,f (x, y),ϕ(x, y) 分别为 Ω,Γ 上的已知连续函数。
第二类和第三类边界条件可统一表示成)0(0),(>=⎪⎭⎫⎝⎛+∂∂Γ∈a u n u y x α (4) 其中 n 为边界 Γ 的外法线方向。
当α = 0 时为第二类边界条件,α ≠ 0 时为第三类边界条件。
在研究热传导过程,气体扩散现象及电磁场的传播等随时间变化的非定常物理问题时,常常会遇到抛物型方程。
其最简单的形式为一维热传导方程)0(022>=∂∂-∂∂a xua t u (5) 方程(5)可以有两种不同类型的定解问题: 初值问题(也称为 Cauchy 问题)⎪⎩⎪⎨⎧+∞<<∞-=+∞<<-∞>=∂∂-∂∂x x x u x t x ua tu )()0,(,0022ϕ (6) 初边值问题⎪⎪⎪⎩⎪⎪⎪⎨⎧<<===<<<<=∂∂-∂∂lx t g t l u t g t u x x u l x T t x ua t u 0),(),(),(),0()()0,(0,002122ϕ (7) 其中ϕ)(),(),(21x g x g x ϕ为已知函数,且满足连接条件 )0()(),0()0(21g l g ==ϕϕ问题(7)中的边界条件)(),(),(),0(21t g t l u t g t u ==称为第一类界条件。
一维抛物线偏微分方程数值解法(4)上一篇参看一维抛物线偏微分方程数值解法(3)(附图及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^4),无条件差分稳定;一:用追赶法解线性方程组(还可以用迭代法解)Matlab程序为:function [u p e x t]=JCHGS(h1,h2,m,n)%紧差分格式解一维抛物线型偏微分方程%此程序用的是追赶法解线性方程组%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)=5/6+r;c(1)=1/12-r/2;d(1)=(r/2-1/12)*u(i+1,1)+... (1/12+r/2)*u(i,1)+(5/6-r)*u(i,2)+(1/12+r/2)*u(i,3)+...h2/12*(f(i,1)+10*f(i,2)+f(i,3));for(k=2:m-2)a(k)=1/12-r/2;b(k)=5/6+r;c(k)=1/12-r/2;d(k)=h2/12*(f(i,k)+...10*f(i,k+1)+f(i,k+2))+(1/12+r/2)*(u(i,k)+u(i,k+2))+(5/6-r)...*u(i,k+1);%输入部分系数矩阵,为0的矩阵元素不输入%一定要注意输入元素的正确性enda(m-1)=1/12-r/2;b(m-1)=5/6+r;d(m-1)=(1/12+r/2)*(u(i,m-1)+u(i,m+1) )+...(5/6-r)*u(i,m)+(r/2-1/12)*u(i+1,m+1)+ ...h2/12*(f(i,m-1)+10*f(i,m)+f(i,m+1));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]=JCHGS(0.1,0.005,10,200); surf(x,t,e)>> title('误差');运行约43秒;[u p e x t]=JCHGS(0.1,0.01,10,100);surf(x,t,e) 20多秒;[u p e x t]=JCHGS(0.2,0.04,5,25);surf(x,t,e) 3秒;此方法精度很高;二:g-s迭代法求解线性方程组Matlab程序function [u e p x t k]=JCFGS1(h1,h2,m,n,kmax,ep) % 解抛物线型一维方程格式(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=((1/12+r/2)*(u(i,j-1)+u(i,j+1))+(5/6-r)*u(i,j)+...h2/12*(f(i,j-1)+10*f(i,j)+f(i,j+1))+(r/2-1/12)*(u(i+1,...j-1)+u(i+1,j+1)))/(5/6+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]=JCFGS1(0.1,0.005,10,200,100000,1e-12);k=67;运行速度1秒左右;surf(x,t,e)[u e p x t k]=JCFGS1(0.01,0.001,100,1000,1000000,1e-12);k=5780;surf(x,t,e)。
引言偏微分方程定解问题有着广泛的应用背景。
人们用偏微分方程来描述、解释或者预见各种自然现象,并用于科学和工程技术的各个领域fll。
然而,对于广大应用工作者来说,从偏微分方程模型出发,使用有限元法或有限差分法求解都要耗费很大的工作量,才能得到数值解。
现在,MATLAB PDEToolbox已实现对于空间二维问题高速、准确的求解过程。
偏微分方程如果一个微分方程中出现的未知函数只含一个自变量,这个方程叫做常微分方程,也简称微分方程;如果一个微分方程中出现多元函数的偏导数,或者说如果未知函数和几个变量有关,而且方程中出现未知函数对几个变量的导数,那么这种微分方程就是偏微分方程。
常用的方法有变分法和有限差分法。
变分法是把定解问题转化成变分问题,再求变分问题的近似解;有限差分法是把定解问题转化成代数方程,然后用计算机进行计算;还有一种更有意义的模拟法,它用另一个物理的问题实验研究来代替所研究某个物理问题的定解。
虽然物理现象本质不同,但是抽象地表示在数学上是同一个定解问题,如研究某个不规则形状的物体里的稳定温度分布问题,由于求解比较困难,可作相应的静电场或稳恒电流场实验研究,测定场中各处的电势,从而也解决了所研究的稳定温度场中的温度分布问题。
随着物理科学所研究的现象在广度和深度两方面的扩展,偏微分方程的应用范围更广泛。
从数学自身的角度看,偏微分方程的求解促使数学在函数论、变分法、级数展开、常微分方程、代数、微分几何等各方面进行发展。
从这个角度说,偏微分方程变成了数学的中心。
一、MATLAB方法简介及应用1.1 MATLAB简介MATLAB是美国MathWorks公司出品的商业数学软件,用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境,主要包括MATLAB和Simulink两大部分。
1.2 Matlab主要功能数值分析数值和符号计算工程与科学绘图控制系统的设计与仿真数字图像处理数字信号处理通讯系统设计与仿真财务与金融工程1.3 优势特点1) 高效的数值计算及符号计算功能,能使用户从繁杂的数学运算分析中解脱出来;2) 具有完备的图形处理功能,实现计算结果和编程的可视化;3) 友好的用户界面及接近数学表达式的自然化语言,使学者易于学习和掌握;4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ,为用户提供了大量方便实用的处理工具。
南京理工大学课程考核论文课程名称:高等数值分析论文题目:有限差分法求解偏微分方程*名:**学号:************成绩:有限差分法求解偏微分方程一、主要内容1.有限差分法求解偏微分方程,偏微分方程如一般形式的一维抛物线型方程:22(,)()u uf x t t xαα∂∂-=∂∂其中为常数具体求解的偏微分方程如下:22001(,0)sin()(0,)(1,)00u u x t x u x x u t u t t π⎧∂∂-=≤≤⎪∂∂⎪⎪⎪=⎨⎪⎪==≥⎪⎪⎩2.推导五种差分格式、截断误差并分析其稳定性;3.编写MATLAB 程序实现五种差分格式对偏微分方程的求解及误差分析;4.结论及完成本次实验报告的感想。
二、推导几种差分格式的过程:有限差分法(finite-difference methods )是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。
有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似,积分用积分和来近似,于是原微分方程和定解条件就近似地代之以代数方程组,即有限差分方程组,解此方程组就可以得到原问题在离散点上的近似解。
推导差分方程的过程中需要用到的泰勒展开公式如下:()2100000000()()()()()()()......()(())1!2!!n n n f x f x f x f x f x x x x x x x o x x n +'''=+-+-++-+- (2-1)求解区域的网格划分步长参数如下:11k k k kt t x x h τ++-=⎧⎨-=⎩ (2-2) 2.1 古典显格式2.1.1 古典显格式的推导由泰勒展开公式将(,)u x t 对时间展开得 2,(,)(,)()()(())i i k i k k k uu x t u x t t t o t t t∂=+-+-∂ (2-3) 当1k t t +=时有21,112,(,)(,)()()(())(,)()()i k i k i k k k k k i k i k uu x t u x t t t o t t tuu x t o tττ+++∂=+-+-∂∂=+⋅+∂ (2-4)得到对时间的一阶偏导数1,(,)(,)()=()i k i k i k u x t u x t uo t ττ+-∂+∂ (2-5) 由泰勒展开公式将(,)u x t 对位置展开得223,,21(,)(,)()()()()(())2!k i k i k i i k i i u uu x t u x t x x x x o x x x x∂∂=+-+-+-∂∂ (2-6)当11i i x x x x +-==和时,代入式(2-6)得2231,1,1122231,1,1121(,)(,)()()()()(())2!1(,)(,)()()()()(())2!i k i k i k i i i k i i i i i k i k i k i i i k i i i iu uu x t u x t x x x x o x x x x u u u x t u x t x x x x o x x x x ++++----⎧∂∂=+-+-+-⎪⎪∂∂⎨∂∂⎪=+-+-+-⎪∂∂⎩(2-7) 因为1k k x x h +-=,代入上式得2231,,22231,,21(,)(,)()()()2!1(,)(,)()()()2!i k i k i k i k i k i k i k i ku uu x t u x t h h o h x xu u u x t u x t h h o h x x +-⎧∂∂=+⋅+⋅+⎪⎪∂∂⎨∂∂⎪=-⋅+⋅+⎪∂∂⎩ (2-8) 得到对位置的二阶偏导数2211,22(,)2(,)(,)()()i k i k i k i k u x t u x t u x t u o h x h+--+∂=+∂ (2-9) 将式(2-5)、(2-9)代入一般形式的抛物线型偏微分方程得21112(,)(,)(,)2(,)(,)(,)()i k i k i k i k i k i k u x t u x t u x t u x t u x t f x t o h h αττ++---+⎡⎤-=++⎢⎥⎣⎦(2-10)为了方便我们可以将式(2-10)写成11122k kk k k k i i i i i i u u u u u f h ατ++-⎡⎤--+-=⎢⎥⎣⎦(2-11) ()11122k k kk k k i i i i i i u u uu u f h τατ++----+= (2-12)最后得到古典显格式的差分格式为()111(12)k k k k k i i i i i u ra u r u u f ατ++-=-+++ (2-13)2r h τ=其中,古典显格式的差分格式的截断误差是2()o h τ+。
《偏微分方程数值解》上机报告实验内容 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 π⎧∂∂=+<<>⎪∂∂⎪⎨==>⎪⎪=+−⎩x 方向 0.1h =, t 方向0.01τ=.在 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.25Current 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时. “向后差分格式”:注:曲线表示精确解,"o"表示数值解(t=0.25时. “六点差分格式” :注:曲线表示精确解,"O"表示数值解(t=0.25时.(三结果分析通过(一 , (二 ,我们检验了三种方法都能很好的求解此一维热传导方程,其中明显能发现“六点对称格式”的误差更小。
第四讲Matlab求解微分方程(组)理论介绍:Matlab求解微分方程(组)命令求解实例:Matlab求解微分方程(组)实例实际应用问题通过数学建模所归纳得到的方程,绝大多数都是微分方程,真正能得到代数方程的机会很少.另一方面,能够求解的微分方程也是十分有限的, 特别是高阶方程和偏微分方程(组).这就要求我们必须研究微分方程(组)的解法:解析解法和数值解法.一.相关函数、命令及简介1.在Matlab中,用大写字母D表示导数,Dy表示y关于自变量的一阶导数, D2y表示y关于自变量的二阶导数,依此类推.函数dsolve用来解决常微分方程(组)的求解问题,调用格式为:X=dsolve(<eqnl,,,eqn2函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解.注意,系统缺省的自变量为t2.函数dsolve求解的是常微分方程的精确解法,也称为常微分方程的符号解. 但是,有大量的常微分方程虽然从理论上讲,其解是存在的,但我们却无法求出其解析解,此时,我们需要寻求方程的数值解,在求常微分方程数值解方面,MATLAB具有丰富的函数,我们将其统称为solver,其一般格式为:[T,Y]=solver(odefun,tspan,yO)说明:(1 )solver 为命令ode45、ode23、odel 13、odel5s、ode23s、ode23t、ode23tb、odel5i 之一.(2)odefun是显示微分方程),=f (t,y)在积分区间tspan =[心心]上从心到“用初始条件儿求解.(3)如果要获得微分方程问题在其他指定时间点bG©…心上的解,则令(span = 『“,•••『/■](要单调的).(4)因为没有一种算法可以有效的解决所有的ODE问题,为此,Matlab提供T多种求解器solver,对于不同的ODE问题,采用不同的solver.程(组)的初值问题的解的Matlab常用程序,其中:ode23采用龙格-库塔2阶算法,用3阶公式作误差估计来调节步长,具有低等的精度.。
二阶偏微分方程的 Matlab有限元法求解摘要:本文基于偏微分方程有限元法求解原理,运用Matlab中的偏微分方程工具箱(PDE Toolbox)对三类典型的二阶偏微分方程:椭圆型方程、双曲线型方程和抛物线型方程算例进行求解,为求解偏微分方程的提供参考。
关键词:偏微分方程,有限元,Matlab偏微分方程工具箱0引言偏微分方程定解问题是描述许多自然现象或工程问题的最重要的数学模型,应用非常广泛[1]。
解析法只能求解非常简单的偏微分方程,远远不能满足科学研究和工程实际的需要。
随着计算机技术和科学计算的迅速发展,数值解法成为求解偏微分方程的重要工具[2-3]。
数值解法将连续问题离散化,最后将偏微分方程化成线性代数方程组。
根据离散化方法不同,偏微分方程数值解法主要有差分法和有限元法。
有限元法是分片定义试函数与变分原理相结合的产物。
它能适应各种形状的区域,且通用性强,现已成为求解偏微分方程定解问题的一种有效数值方法[4]。
本文首先简述了偏微分方程有限元法原理,然后,对Matlab中的偏微分方程工具箱(Partial Differential Equations Toolbox)的功能和求解思路进行了阐述[5-6],最后,给出了用PDE Toolbox求解椭圆方程、、双曲线方程和抛物线方程的计算实例。
1偏微分方程有限元法原理偏微分方程有限元法的基本思想是将实际上连续的整个求解域进行离散化处理,即用一些假想的面或线将求解域分割为一系列的单元,各个单元之间仅在有限个节点处相互连接。
取未知函数的节点值作为基本未知量,在每个单元上选取一个近似的插值函数表示单元中场函数的分布规律。
利用变分原理来获得单元的刚度方程,然后按一定的规则把所有单元的刚度方程组集合起来,经适当的边界条件处理,便得到整个系统的总体方程组。
这样,偏微分方程便转化为一组常微分方程。
最后,求解总体方程组,得到节点值和用插值函数确定整个求解域上的场函数。
matlab解偏微分方程Matlab是一种非常强大的数学计算工具,它可以用于解决各种数学问题。
在本文中,我们将学习如何使用Matlab解偏微分方程。
偏微分方程是一类包含未知函数的偏导数的方程。
通常,解偏微分方程是困难的,需要使用复杂的数学方法。
然而,Matlab可以大大简化这个过程。
在Matlab中,我们可以使用pdepe函数来解偏微分方程。
pdepe函数采用一个偏微分方程的系统,并返回一个包含解的向量的矩阵。
下面是一个解二维扩散方程的示例程序:%定义二维扩散方程 function [c,f,s] = diffusionpde(x,t,u,DuDx)c = 1; %系数f = DuDx; %带有时间和空间导数的项s = 0; %不带导数的项end%定义边界条件(例)function [pl,ql,pr,qr] =diffusionbc(xl,ul,xr,ur,t)pl = 0; ql = 1; %左边界(u=0)pr = 0; qr = 1; %右边界(u=0)end%定义初始条件(例)function u0 = diffusionic(x)u0 = sin(pi*x); %sin(pi*x)是初始条件方程end%主程序x = linspace(0,1,50); %空间网格t = linspace(0,1,10); %时间网格sol =pdepe(0,@diffusionpde,@diffusionic,@diffusionbc,x,t );u = sol(:,:,1); %提取第一个解%绘制解surfc(x,t,u)xlabel('位置')ylabel('时间')title('二维扩散方程的解')从上述程序中,我们可以看到pdepe的使用方法。
在主程序中,我们选择了空间和时间网格,然后定义了偏微分方程、初始条件和边界条件的函数。
最后,我们调用pdepe函数,并将解存储在变量sol中。
MATLAB中的偏微分方程数值解法偏微分方程(Partial Differential Equations,PDEs)是数学中的重要概念,广泛应用于物理学、工程学、经济学等领域。
解决偏微分方程的精确解往往非常困难,因此数值方法成为求解这类问题的有效途径。
而在MATLAB中,有丰富的数值解法可供选择。
本文将介绍MATLAB中几种常见的偏微分方程数值解法,并通过具体案例加深对其应用的理解。
一、有限差分法(Finite Difference Method)有限差分法是最为经典和常用的偏微分方程数值解法之一。
它将偏微分方程的导数转化为差分方程,通过离散化空间和时间上的变量,将连续问题转化为离散问题。
在MATLAB中,使用有限差分法可以比较容易地实现对偏微分方程的数值求解。
例如,考虑一维热传导方程(Heat Equation):∂u/∂t = k * ∂²u/∂x²其中,u为温度分布随时间和空间的变化,k为热传导系数。
假设初始条件为一段长度为L的棒子上的温度分布,边界条件可以是固定温度、热交换等。
有限差分法可以将空间离散化为N个节点,时间离散化为M个时刻。
我们可以使用中心差分近似来计算二阶空间导数,从而得到以下差分方程:u(i,j+1) = u(i,j) + Δt * (k * (u(i+1,j) - 2 * u(i,j) + u(i-1,j))/Δx²)其中,i表示空间节点,j表示时间步。
Δt和Δx分别为时间和空间步长。
通过逐步迭代更新节点的温度值,我们可以得到整个时间范围内的温度分布。
而MATLAB提供的矩阵计算功能,可以大大简化有限差分法的实现过程。
二、有限元法(Finite Element Method)有限元法是另一种常用的偏微分方程数值解法,特点是适用于复杂的几何形状和边界条件。
它将求解区域离散化为多个小单元,通过构建并求解代数方程组来逼近连续问题。
在MATLAB中,我们可以使用Partial Differential Equation Toolbox提供的函数进行有限元法求解。