有限元解二维变系数抛物方程
- 格式:doc
- 大小:22.00 KB
- 文档页数:3
抛物型方程的Galerkin有限元方法一、引言抛物型方程是一类常见的偏微分方程,具有广泛的应用。
在数值解中,Galerkin有限元方法是一种常用且有效的方法。
本文将介绍抛物型方程的基本概念,并详细讲解Galerkin有限元方法在求解抛物型方程中的应用。
二、抛物型方程的基本概念抛物型方程是指具有二阶时间导数和二阶空间导数的偏微分方程。
一般形式为:∂u−Δu=f∂t其中,u为未知函数,t为时间变量,Δ为Laplace算子,f为给定的函数。
抛物型方程的一个重要特点是初始条件和边界条件对解的影响非常大。
合适的初始条件和边界条件能够唯一确定方程的解。
三、Galerkin有限元方法Galerkin有限元方法是一种利用函数空间进行近似的数值计算方法。
它基于以下思想:将问题的解表示为函数空间中的一个函数,通过求解一组代数方程组来近似求解原始方程。
1. 函数空间的选择在应用Galerkin有限元方法求解抛物型方程时,需要选择合适的函数空间。
常用的函数空间有有限维函数空间和无限维函数空间。
具体的选择需要根据问题的特点和计算的要求来确定。
2. 弱形式的推导对于抛物型方程,我们可以将其转化为弱形式。
弱形式是通过将方程两边乘以一个测试函数,并进行积分得到的。
这样可以减小对解的要求,并使得问题更容易求解。
3. 数值离散和代数方程的建立接下来,需要对时间和空间进行离散。
通常使用网格来进行离散,将时间和空间分割为有限个小区域。
然后,通过选择适当的基函数,在每个小区域上近似原方程的解。
最终得到一组代数方程组。
求解代数方程组是Galerkin有限元方法的最后一步。
可以使用常用的数值方法,如迭代法、直接法等,来求解代数方程组。
根据计算要求和问题特点,选择合适的求解方法。
四、应用案例以一维热传导方程为例,展示Galerkin有限元方法在求解抛物型方程中的应用。
热传导方程是一个典型的抛物型方程,描述了物体内部的温度分布随时间变化的规律。
《两类方程的时空混合有限元格式》篇一一、引言在数值分析和计算物理领域,有限元方法是一种重要的数值技术,广泛应用于解决各类复杂的偏微分方程问题。
近年来,随着对时空混合问题研究的深入,混合有限元格式成为了一种具有广泛应用前景的数值求解方法。
本文将介绍两种不同类型的方程,即抛物型方程和双曲型方程的时空混合有限元格式,并探讨其在实际问题中的应用。
二、抛物型方程的时空混合有限元格式抛物型方程是描述物质在时间上的扩散、热传导等过程的数学模型。
对于这类问题,我们采用时空混合有限元方法进行求解。
首先,我们将时间域和空间域进行离散化处理,然后通过构建适当的基函数,将原问题转化为一个线性系统。
接着,利用高斯消元法或迭代法求解该线性系统,得到问题的解。
在抛物型方程的时空混合有限元格式中,我们采用了等参元和线性插值技术,使得求解过程更加稳定和高效。
此外,我们还采用了自适应网格技术,根据问题的特点自动调整网格的疏密程度,进一步提高求解精度。
三、双曲型方程的时空混合有限元格式双曲型方程是描述物质在时间和空间上传播、振动等过程的数学模型。
与抛物型方程类似,我们同样采用时空混合有限元方法进行求解。
然而,由于双曲型方程具有更复杂的动力学特性,我们需要采用更加精细的离散化方法和基函数来描述问题的解。
在双曲型方程的时空混合有限元格式中,我们采用了高阶基函数和精细的网格划分来提高求解精度。
同时,我们还采用了稳定化技术来处理数值计算中的不稳定性问题。
此外,我们还探讨了多尺度问题的处理方法,通过引入多尺度基函数来提高求解效率。
四、实际应用我们通过两个实际问题的求解过程来展示两类方程的时空混合有限元格式的应用。
首先是一个二维的热传导问题,我们利用抛物型方程的时空混合有限元格式求解了该问题,得到了良好的数值结果。
其次是一个弹丸在空气中传播的问题,我们利用双曲型方程的时空混合有限元格式进行了模拟和分析,得到了弹丸传播过程中的速度、压力等物理量的变化情况。
分类号:O241.82本科生毕业论文(设计)题目:一类抛物型方程的计算方法作者单位数学与信息科学学院作者姓名专业班级2011级数学与应用数学创新2班指导教师论文完成时间二〇一五年四月一类抛物型方程的数值计算方法(数学与信息科学学院数学与应用数学专业2011级创新2班)指导教师摘要: 抛物型方程数值求解常用方法有差分方法、有限元方法等。
差分方法是一种对方程直接进行离散化后得到的差分计算格式,有限元方法是基于抛物型方程的变分形式给出的数值计算格式。
本文首先给出抛物型方程的差分计算方法,并分析了相应差分格式的收敛性、稳定性等基本理论问题.然后,给出抛物型方程的有限元计算方法及理论分析。
关键词:差分方法,有限元方法,收敛性,稳定性Numerical computation methods for a parabolic equationYan qian(Class 2, Grade 2011,College of Mathematics and Information Science)Advisor: Nie huaAbstract:The common methods to solve parabolic equations include differential method,finite element method etc。
The main idea of differential method is to construct differential schemes by discretizing differential equations directly. Finite element scheme is based on the variational method of parabolic equations。
In this article, we give some differential schemes for a parabolic equation and analyze their convergence and stability. Moreover,the finite element method and the corresponding theoretical analysis for parabolic equation are established.Key words:differential method,finite element method, convergence,stability1 绪 论1。
ADI隐式交替法三种解法及误差分析(一般的教材上只说第一种)理论部分参看孙志忠:偏微分方程数值解法注意:1.最好不要直接看程序,中间很多公式很烦人的(一定要小心),我写了两天,终于写对了。
2.中间:例如r*(u(i-1,m1,k)+u(i+1,m1,k))形式写成分形式:r*u(i-1,m1,k)+r*u(i+1,m1,k)后面会出错,我也不是很清楚为什么,可能由于舍入误差,或者大数吃掉小数的影响。
3.下面有三个程序4.具体理论看书,先仔细看书(孙志忠:偏微分方程数值解法)或者网上搜一些理论。
Matlab程序:1.function [u u0 p e x y t]=ADI1(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(P-R交替隐式,截断)%此程序用的是追赶法解线性方程组%h1为空间步长,h2为时间步长%m1,m2分别为x方向,y方向网格数,n为时间网格数%p为精确解,u为数值解,e为误差%定义u0(i,j,k)=u(i,j,k+1/2),因为矩阵中,i,j,k必须全为整数x=(0:m1)*h1+0;%定义x0,y0,t0是为了f(x,t)~=0的情况%y=(0:m2)*h1+0;t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;for k=1:n+1for i=1:m2+1for j=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));endendendfor i=1:m2+1for j=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endendfor k=1:n+1for i=1:m2+1u(i,[1 m1+1],k)=[exp(0.5*y(i)-t(k)) exp(0.5*(1+y(i))-t(k))]; u0(i,[1 m1+1],k)=[exp(0.5*y(i)-t0(k)) exp(0.5*(1+y(i))-t0(k))] ;endendfor k=1:n+1for j=1:m1+1u([1 m2+1],j,k)=[exp(0.5*x(j)-t(k)) exp(0.5*(1+x(j))-t(k))]; u0([1 m2+1],j,k)=[exp(0.5*x(j)-t0(k)) exp(0.5*(1+x(j))-t0(k))];endendr=h2/(h1*h1);r1=2*(1-r);r2=2*(1+r);for k=1:n %外循环,先固定每一时间层,每一时间层上解一线性方程组% for i=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=r2*ones(1,m1-1);d(1)=r*u0(i,1,k)+r*(u(i-1,2,k)+u(i+1,2,k))+r1*u(i,2,k)+...h2*f(i,2,k);for l=2:m1-2d(l)=r*(u(i-1,l+1,k)+u(i+1,l+1,k))+r1*u(i,l+1,k)+...h2*f(i,l+1,k);%输入部分系数矩阵,为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r*(u(i-1,m1,k)+u(i+1,m1,k))...+r1*u(i,m1,k)+h2*f(i,m1,k);for l=1:m1-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu0(i,m1,k)=d(m1-1)/b(m1-1); %回代过程%for l=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendfor j=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=r2*ones(1,m2-1);d(1)=r*u(1,j,k+1)+r*(u0(2,j-1,k)+u0(2,j+1,k))+r1*u0(2,j,k)+...h2*f(2,j,k);for l=2:m2-2d(l)=r*(u0(l+1,j-1,k)+u0(l+1,j+1,k))+r1*u0(l+1,j,k)+... h2*f(l+1,j,k);%输入部分系数矩阵,为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=r*u(m2+1,j,k+1)+r*(u0(m2,j-1,k)+u0(m2,j+1,k))...+r1*u0(m2,j,k)+h2*f(m2,j,k);for l=1:m2-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%for l=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfor k=1:n+1for i=1:m2+1for j=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k)); %p为精确解e(i,j,k)=abs(u(i,j,k)-p(i,j,k)); %e为误差endendend2.function [u p e x y t]=ADI2(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(D'Yakonov交替方向隐格式)%此程序用的是追赶法解线性方程组%h1为空间步长,h2为时间步长%m1,m2分别为x方向,y方向网格数,n为时间网格数%p为精确解,u为数值解,e为误差%定义u0(i,j,k)=u'(i,j,k)(引入的过渡层),因为矩阵中,i,j,k必须全为整数x=(0:m1)*h1+0;y=(0:m2)*h1+0;t=(0:n)*h2+0;t0=(0:n)*h2+1/2*h2;%定义t0是为了f(x,y,t)~=0的情况% for k=1:n+1for i=1:m2+1for j=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));%编程时-t0(k)写成了+t0(k),导致错误;endendend%初始条件for i=1:m2+1for j=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endend%边界条件for k=1:n+1for i=1:m2+1u(i,[1 m1+1],k)=[exp(0.5*y(i)-t(k)) exp(0.5*(1+y(i))-t(k))];endendr=h2/(h1*h1);r4=1+r;r5=r/2;for k=1:nfor i=2:m2u0(i,[1 m1+1],k)=r4*u(i,[1 m1+1],k+1)-r5*(u(i-1,[1 m1+1],... k+1)+u(i+1,[1 m1+1],k+1));endendfor k=1:n+1for j=1:m1+1u([1 m2+1],j,k)=[exp(0.5*x(j)-t(k)) exp(0.5*(1+x(j))-t(k))];endendr1=r-r*r;r2=2*(r-1)*(r-1);r3=r*r/2;for k=1:n %外循环,先固定每一时间层,每一时间层上解一线性方程组% for i=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=2*r4*ones(1,m1-1);d(1)=r*u0(i,1,k)+r1*(u(i-1,2,k)+u(i,1,k)+u(i+1,2,k)+...u(i,3,k))+r2*u(i,2,k)+r3*(u(i-1,1,k)+...u(i+1,1,k)+u(i-1,3,k)+u(i+1,3,k))+2*h2*f(i,2,k);for l=2:m1-2d(l)=r1*(u(i-1,l+1,k)+u(i,l,k)+u(i+1,l+1,k)+...u(i,l+2,k))+r2*u(i,l+1,k)+r3*(u(i-1,l,k)+...u(i+1,l,k)+u(i-1,l+2,k)+u(i+1,l+2,k))+2*h2*f(i,l+1,k);%输入部分系数矩阵,为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r1*(u(i-1,m1,k)+u(i,m1-1,k)+...u(i+1,m1,k)+u(i,m1+1,k))+r2*u(i,m1,k)+...r3*(u(i-1,m1-1,k)+...u(i+1,m1-1,k)+u(i-1,m1+1,k)+u(i+1,m1+1,k))+2*h2*f(i,m1,k);for l=1:m1-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);end%回代过程%u0(i,m1,k)=d(m1-1)/b(m1-1);for l=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendfor j=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=2*r4*ones(1,m2-1);d(1)=r*u(1,j,k+1)+2*u0(2,j,k);for l=2:m2-2d(l)=2*u0(l+1,j,k);%输入部分系数矩阵,为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=2*u0(m2,j,k)+r*u(m2+1,j,k+1);for l=1:m2-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%for l=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfor k=1:n+1for i=1:m2+1for j=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k)); %p为精确解e(i,j,k)=abs(u(i,j,k)-p(i,j,k)); %e为误差endendend3.function [u u0 p e x y t]=ADI5(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(P-R交替隐式,未截断)%此程序用的是追赶法解线性方程组%h1为空间步长,h2为时间步长%m1,m2分别为x方向,y方向网格数,n为时间网格数%p为精确解,u为数值解,e为误差%定义u0(i,j,k)=u(i,j,k+1/2),因为矩阵中,i,j,k必须全为整数x=(0:m1)*h1+0;%定义x0,y0,t0是为了f(x,t)~=0的情况%y=(0:m2)*h1+0;t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;for k=1:n+1for i=1:m2+1for j=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));endendendfor i=1:m2+1for j=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endendfor k=1:n+1for i=1:m2+1u(i,[1 m1+1],k)=[exp(0.5*y(i)-t(k)) exp(0.5*(1+y(i))-t(k))]; u1(i,[1 m1+1],k)=[exp(0.5*y(i)-t0(k)) exp(0.5*(1+y(i))-t0(k))] ;endendr=h2/(h1*h1);r1=2*(1-r);r2=r/4;r3=2*(1+r);for k=1:nfor i=2:m2u0(i,[1 m1+1],k)=u1(i,[1 m1+1],k)-r2*(u(i-1,[1 m1+1],k+1)-...2*u(i,[1 m1+1],k+1)+u(i+1,[1 m1+1],k+1)-u(i-1,[1 m1+1],k)+...2*u(i,[1 m1+1],k)-u(i+1,[1 m1+1],k));endendfor k=1:n+1for j=1:m1+1u([1 m2+1],j,k)=[exp(0.5*x(j)-t(k)) exp(0.5*(1+x(j))-t(k))];endendfor k=1:n %外循环,先固定每一时间层,每一时间层上解一线性方程组%for i=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=r3*ones(1,m1-1);d(1)=r*u0(i,1,k)+r*(u(i-1,2,k)+u(i+1,2,k))+r1*u(i,2,k)+... h2*f(i,2,k);for l=2:m1-2d(l)=r*(u(i-1,l+1,k)+u(i+1,l+1,k))+r1*u(i,l+1,k)+...h2*f(i,l+1,k);%输入部分系数矩阵,为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r*(u(i-1,m1,k)+u(i+1,m1,k))...+r1*u(i,m1,k)+h2*f(i,m1,k);for l=1:m1-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu0(i,m1,k)=d(m1-1)/b(m1-1); %回代过程%for l=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendfor j=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=r3*ones(1,m2-1);d(1)=r*u(1,j,k+1)+r*(u0(2,j-1,k)+u0(2,j+1,k))+r1*u0(2,j,k)+...h2*f(2,j,k);for l=2:m2-2d(l)=r*(u0(l+1,j-1,k)+u0(l+1,j+1,k))+r1*u0(l+1,j,k)+... h2*f(l+1,j,k);%输入部分系数矩阵,为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=r*u(m2+1,j,k+1)+r*(u0(m2,j-1,k)+u0(m2,j+1,k))...+r1*u0(m2,j,k)+h2*f(m2,j,k);for l=1:m2-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%for l=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfor k=1:n+1for i=1:m2+1for j=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k)); %p为精确解 e(i,j,k)=abs(u(i,j,k)-p(i,j,k)); %e为误差endendend[up e x y t]=ADI2(0.01,0.001,100,100,1000);surf(x,y,e(:,:,1001)) t=1的误差曲面下面是三种方法的误差比较:1.[u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式,截断)截断中间过渡层用u(i,j,k+1/2)代替)(t=1时的误差)2.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式,未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)3.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式) surf(x,y,e(:,:,11))(表示t=1时的误差)下面是相关数据:1: [u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式,截断)截断中间过渡层用u(i,j,k+1/2)代替)e(:,:,11) =Columns 1 through 60 0 0 0 0 00 0.00040947 0.00025182 0.00019077 0.00017112 0.000176040 0.00057359 0.00042971 0.00035402 0.00032565 0.000336280 0.00066236 0.00054689 0.00047408 0.00044596 0.000462670 0.00072152 0.00062001 0.00055081 0.00052442 0.000545530 0.00076164 0.0006576 0.00058522 0.00055732 0.000579840 0.00078336 0.00065993 0.00057557 0.00054161 0.000562090 0.00078161 0.00061872 0.00051646 0.00047429 0.000489640 0.00073621 0.0005148 0.00039979 0.00035439 0.000363130 0.00056964 0.00031688 0.00022051 0.0001884 0.000191920 0 0 0 0 02.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式,未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)e(:,:,11) =Columns 1 through 60 0 0 0 0 00 0.00027006 0.00016305 0.00012104 0.0001071 0.000109950 0.00037754 0.00027817 0.0002253 0.00020483 0.000211160 0.00043539 0.00035386 0.00030207 0.00028124 0.00029140 0.00047398 0.00040104 0.00035113 0.00033111 0.000344050 0.0005003 0.00042535 0.00037309 0.0003519 0.000365710 0.00051479 0.00042699 0.00036681 0.00034164 0.00035410 0.00051415 0.00040056 0.00032887 0.0002985 0.000307640 0.00048504 0.0003335 0.00025411 0.0002221 0.000227060 0.00037609 0.00020532 0.00013956 0.00011718 0.000119020 0 0 0 0 03.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式)e(:,:,11) =Columns 1 through 60 0 0 0 0 00 8.6469e-006 1.4412e-005 1.8364e-005 2.091e-005 2.2174e-0050 1.4412e-005 2.4777e-005 3.2047e-005 3.6716e-005 3.8961e-0050 1.8364e-005 3.2047e-005 4.1789e-005 4.8054e-005 5.1008e-0050 2.091e-005 3.6716e-005 4.8054e-005 5.5353e-005 5.8764e-0050 2.2174e-005 3.8961e-005 5.1008e-005 5.8764e-005 6.2389e-0050 2.2118e-005 3.8698e-005 5.0523e-005 5.8126e-005 6.171e-0050 2.055e-005 3.5581e-005 4.6157e-005 5.2942e-005 5.6197e-0050 1.707e-005 2.8951e-005 3.7128e-005 4.2365e-005 4.4952e-0050 1.0851e-005 1.7698e-005 2.2265e-005 2.5203e-005 2.672e-0050 0 0 0 0 01.[u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式,截断)截断中间过渡层用u(i,j,k+1/2)代替)Columns 7 through 110 0 0 0 00.00020348 0.00026228 0.00038338 0.00066008 00.00038607 0.00048321 0.00064717 0.00091668 00.00052635 0.00064203 0.00081637 0.0010517 00.0006174 0.00074272 0.00092111 0.0011417 00.00065651 0.00078964 0.00097724 0.0012051 00.00064051 0.00078116 0.00098594 0.0012433 00.00056474 0.00070822 0.00093332 0.0012478 00.00042547 0.00055526 0.00078616 0.0011844 00.00022735 0.00030946 0.00049004 0.00092402 00 0 0 0 02.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式,未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)Columns 7 through 110 0 0 0 00.00012826 0.00016798 0.00024986 0.00043637 00.00024444 0.00031023 0.00042173 0.00060513 00.00033401 0.00041257 0.00053179 0.00069358 00.00039216 0.00047742 0.00059986 0.00075263 00.00041704 0.00050761 0.00063642 0.00079439 00.00040657 0.0005021 0.00064226 0.00081984 00.00035784 0.000455 0.00060828 0.00082334 00.00026866 0.00035628 0.00051263 0.0007824 00.00014262 0.00019789 0.00031956 0.0006113 00 0 0 0 0 3.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式) Columns 7 through 110 0 0 0 02.2118e-005 2.055e-005 1.707e-005 1.0851e-005 03.8698e-005 3.5581e-005 2.8951e-005 1.7698e-005 05.0523e-005 4.6157e-005 3.7128e-005 2.2265e-005 05.8126e-005 5.2942e-005 4.2365e-005 2.5203e-005 06.171e-005 5.6197e-005 4.4952e-005 2.672e-005 06.1116e-005 5.5803e-005 4.4817e-005 2.6785e-005 05.5803e-005 5.1239e-005 4.1529e-005 2.5153e-005 04.4817e-005 4.1529e-005 3.42e-005 2.126e-005 02.6785e-005 2.5153e-005 2.126e-005 1.3869e-005 00 0 0 0 0。
二维抛物方程的有限差分法摘要二维抛物方程是一类有广泛应用的偏微分方程,由于大部分抛物方程都难以求得解析解,故考虑采用数值方法求解。
有限差分法是最简单又极为重要的解微分方程的数值方法。
本文介绍了二维抛物方程的有限差分法。
首先,简单介绍了抛物方程的应用背景,解抛物方程的常见数值方法,有限差分法的产生背景和发展应用。
讨论了抛物方程的有限差分法建立的基础,并介绍了有限差分方法的收敛性和稳定性。
其次,介绍了几种常用的差分格式,有古典显式格式、古典隐式格式、Crank-Nicolson隐式格式、Douglas差分格式、加权六点隐式格式、交替方向隐式格式等,重点介绍了古典显式格式和交替方向隐式格式。
进行了格式的推导,分析了格式的收敛性、稳定性。
并以热传导方程为数值算例,运用差分方法求解。
通过数值算例,得出古典显式格式计算起来较简单,但稳定性条件较苛刻;而交替方向隐式格式无条件稳定。
关键词:二维抛物方程;有限差分法;古典显式格式;交替方向隐式格式FINITE DIFFERENCE METHOD FORTWO-DIMENSIONAL PARABOLICEQUATIONAbstractTwo-dimensional parabolic equation is a widely used class of partial differential equations. Because this kind of equation is so complex, we consider numerical methods instead of obtaining analytical solutions. finite difference method is the most simple and extremely important numerical methods for differential equations. The paper introduces the finite difference method for two-dimensional parabolic equation.Firstly, this paper introduces the background and common numerical methods for Parabolic Equation, Background and development of applications. Discusses the basement for the establishment of the finite difference method for parabolic equation And describes the convergence and stability for finite difference method.Secondly, Introduces some of the more common simple differential format,for example, the classical explicit scheme, the classical implicit scheme, Crank-Nicolson implicit scheme, Douglas difference scheme, weighted six implicit scheme and the alternating direction implicit format. The paper focuses on the classical explicit scheme and the alternating direction implicit format. The paper takes discusses the derivation convergence,and stability of the format . The paper takes And the heat conduction equation for the numerical example, using the differential method to solve. Through numerical examples, the classical explicit scheme is relatively simple for calculation, with more stringent stability conditions; and alternating direction implicit scheme is unconditionally stable.Keywords:Two-dimensional Parabolic Equation; Finite-Difference Method; Eclassical Explicit Scheme; Alternating Direction Implicit Scheme目录摘要 (I)Abstract (II)1绪论 (1)1.1课题背景 (1)1.2发展概况 (1)1.2.1抛物型方程的常见数值解法 (1)1.2.2有限差分方法的发展 (2)1.3差分格式建立的基础 (3)1.3.1区域剖分 (3)1.3.2差商代替微商 (3)1.3.3差商代替微商格式的误差分析 (4)1.4本文主要研究容 (5)2显式差分格式 (7)2.1常系数热传导方程的古典显式格式 (7)2.1.1古典显式格式格式的推导 (7)2.1.3古典显式格式的算法步骤 (8)3隐式差分格式 (10)3.1古典隐式格式 (10)3.2 Crank-Nicolson隐式格式 (12)3.3 Douglas差分格式 (13)3.4加权六点隐式格式 (14)3.5交替方向隐式格式 (15)3.5.1 Peaceman-Rachford格式 (15)3.5.2 Rachford-Mitchell格式 (15)3.5.3 Mitchell-Fairweather格式 (15)3.5.4交替方向隐式格式的算法步骤 (16)4实例分析与结果分析 (17)4.1算例 (17)4.1.1已知有精确解的热传导问题 (17)4.1.2未知精确解的热传导问题 (19)4.2结果分析 (20)5稳定性探究与分析 (21)5.1稳定性问题的提出 (21)5.2 几种分析稳定性的方法 (21)5.3 r变化对稳定性的探究 (23)5.3.1 古典显式格式的稳定性 (23)5.3.2 P-R格式格式的稳定性 (24)结语 (26)参考文献 (27)附录P-R格式的C++实现代码 (28)致谢 (30)1绪论1.1课题背景抛物方程是一类特殊的偏微分方程,二维抛物方程的一般形式为u Lu t∂=∂ (1-1) 其中1212((,,))((,,))(,,)(,,)(,,)u u u u u u L a x y t a x y t b x y t b x y t C x y t x x y y x y∂∂∂∂∂∂=++++∂∂∂∂∂∂ 120,0,0a a C >>≥。
python有限元差分求解二维偏微分方程摘要:1.引言2.有限元方法简介3.差分方法简介4.二维偏微分方程的求解方法5.Python 有限元差分求解二维偏微分方程的实现6.总结正文:【引言】本文旨在介绍如何使用Python 有限元差分法求解二维偏微分方程。
有限元方法和差分方法是数值计算领域中常用的方法,它们可以用来解决偏微分方程问题。
在本文中,我们将使用Python 编程语言来实现这些方法,以解决二维偏微分方程。
【有限元方法简介】有限元方法是一种数值计算方法,它将连续的求解区域离散化为有限个单元,通过在每个单元内求解局部问题,最后将各单元的解合并得到原问题的解。
有限元方法广泛应用于固体力学、流体力学、热传导等领域。
【差分方法简介】差分方法是一种数值计算方法,它通过将连续的函数值用离散的点表示,然后通过差分公式将函数在某点的导数表示为有限差分形式。
差分方法主要包括前向差分、后向差分和中心差分等。
差分方法在数值计算中具有重要地位,例如在求解常微分方程的初值问题、边值问题以及偏微分方程等方面都有广泛应用。
【二维偏微分方程的求解方法】二维偏微分方程的求解方法主要包括有限元方法和差分方法。
有限元方法将求解区域离散化为有限个单元,通过在每个单元内求解局部问题,最后将各单元的解合并得到原问题的解。
而差分方法通过将连续的函数值用离散的点表示,然后通过差分公式将函数在某点的导数表示为有限差分形式。
这两种方法在求解二维偏微分方程时,可以相互结合,提高求解精度。
【Python 有限元差分求解二维偏微分方程的实现】Python 作为一门功能强大的编程语言,提供了丰富的数值计算库,如NumPy 和SciPy。
利用这些库,我们可以方便地实现有限元差分法求解二维偏微分方程。
下面是一个简单的示例,使用Python 求解二维Laplace 方程:```pythonimport numpy as npfrom scipy.spatial import griddatadef Laplace(x, y):return np.sin(x) + np.cos(y)def Laplace_dirichlet(x, y, x0, y0):return np.sin(x) + np.cos(y) - np.sin(x0) - np.cos(y0) # 网格划分x = np.linspace(0, 2 * np.pi, 100)y = np.linspace(0, 2 * np.pi, 100)X, Y = np.meshgrid(x, y)# 边界条件x0, y0 = 0, np.pi / 2# 求解Z = griddata(X, Y, Laplace(X, Y), (x, y), method="cubic")Z_dirichlet = griddata(X, Y, Laplace_dirichlet(X, Y, x0, y0), (x, y), method="cubic")# 绘制结果import matplotlib.pyplot as pltplt.figure()plt.contourf(x, y, Z)plt.colorbar()plt.title("二维Laplace 方程的数值解")plt.show()plt.figure()plt.contourf(x, y, Z_dirichlet)plt.colorbar()plt.title("二维Laplace 方程的有限元差分解")plt.show()```【总结】本文介绍了Python 有限元差分求解二维偏微分方程的方法。
抛物型方程的galerkin有限元方法抛物型方程是一类重要的偏微分方程,它在物理、工程、经济等领域中都有广泛的应用。
而galerkin有限元方法是一种常用的数值解法,可以有效地求解抛物型方程。
本文将介绍抛物型方程的galerkin有限元方法。
一、抛物型方程抛物型方程是一类偏微分方程,其一般形式为:$$\frac{\partial u}{\partial t} - \nabla \cdot (a\nabla u) + cu = f $$其中,$u$是未知函数,$a$和$c$是已知函数,$f$是给定函数。
抛物型方程的特点是时间和空间都是连续的,因此需要使用时间和空间上的离散化方法来求解。
二、galerkin有限元方法galerkin有限元方法是一种常用的数值解法,它将偏微分方程的解表示为一组基函数的线性组合,然后通过求解系数来得到解。
具体来说,galerkin有限元方法将偏微分方程的解表示为:$$u_h(x,t) = \sum_{i=1}^N u_i(t) \phi_i(x)$$其中,$u_i(t)$是待求系数,$\phi_i(x)$是一组基函数,$N$是基函数的个数。
将上式代入偏微分方程中,得到:$$\sum_{i=1}^N \left( \frac{\partial u_i}{\partial t} - \nabla \cdot(a\nabla \phi_i) + c\phi_i \right) \phi_j = \int_\Omega f\phi_j $$对于任意的$j=1,2,\cdots,N$,上式都成立。
因此,可以得到一个关于系数$u_i(t)$的线性方程组,通过求解该方程组即可得到解$u_h(x,t)$。
三、抛物型方程的galerkin有限元方法将抛物型方程代入galerkin有限元方法中,得到:\sum_{i=1}^N \left( \frac{\partial u_i}{\partial t} - \nabla \cdot (a\nabla \phi_i) + c\phi_i \right) \phi_j = \int_\Omega f\phi_j $$对于任意的$j=1,2,\cdots,N$,上式都成立。
二维抛物线方程数值解法方法二维抛物线方程是数学中常见的偏微分方程之一,表示一个二维区域中的传热、扩散、波动等现象。
求解二维抛物线方程通常需要使用数值解法,其中ADI(Alternating Direction Implicit)隐式交替法是一种有效的数值解法之一、本文将详细介绍ADI隐式交替法的原理和步骤。
首先,我们来看一下二维抛物线方程的一般形式:∂u/∂t=α(∂²u/∂x²+∂²u/∂y²)+f(x,y,t)其中u(x,y,t)是未知函数,表示二维区域中的一些物理量(如温度),α是扩散系数,f(x,y,t)是源(即外力项)。
ADI隐式交替法的基本思想是将偏微分方程中的时间导数项进行分离,然后通过交替方向的半隐式差分逼近来求解。
具体步骤如下:1.将时间导数项进行分离,得到两个互相独立的方程:∂u/∂t=α(∂²u/∂x²)+α(∂²u/∂y²)+f(x,y,t)等价于:∂u/∂t=α(∂²u/∂x²)+f(x,y,t),考虑到y方向上的导数已经被分离出来。
∂u/∂t=α(∂²u/∂y²)+f(x,y,t),考虑到x方向上的导数已经被分离出来。
2. 对第一个方程应用隐式差分(比如Crank-Nicolson差分格式),得到一个关于未知函数u的线性方程组Ax=b。
3.利用ADI思想,对第二个方程同样进行隐式差分,得到另一个关于未知函数u的线性方程组Ay=c。
4.对线性方程组Ax=b和Ay=c进行求解,得到u的一个时间层的近似解。
5.重复步骤2至4,直到得到整个时间区间内的数值解。
ADI隐式交替法的关键在于递归地求解两个互相独立的线性方程组。
在每个时间步长中,先求解一个方向(比如x方向)上的方程组,然后用该方程组的解代入另一个方向(比如y方向)的方程组中求解。
通过交替进行,可以逼近二维抛物线方程的解。
前言抛物型方程解的估计及其应用1前言数学物理方程主要指从物理学及其它各门自然科学、技术科学中所产生的偏微分方程(有时也包括积分方程、微分方程等),它们反映了有关的未知变量关于时间的导数和关于空间变量的导数之间的制约关系.连续介质力学、电磁学、量子力学等等方面的基本方程属于数学物理方程的范围.它以具有物理背景的偏微分方程(组)作为研究的主要对象.它与其他数学分支及物理、化学等自然科学和工程技术的很多领域都有着广泛的联系,因此,无论在历史上还是在今天的现实生活中,它对于推动数学理论的发展,加强理论与实际的联系,帮助人们认识世界和改造世界都起着重要大的作用.微积分产生以后,人们就开始把力学中的一些问题,归结为偏微分方程进行研究.早在18世纪初,人们已经将弦线振动问题归结为弦振动方程,并探讨了它的解法.随后,人们又陆续了解了流体的运动、弹性体的平衡和振动、热传导、电磁相互作用、原子核和电子的相互作用、化学反应过程等等自然现象的基本规律,把它们写成偏微分方程的形式,并且求出了典型问题的解答,从而能通过实践,验证这些基本规律的正确性,显示了数学物理方程对于认识自然界基本规律的重要性.有了基本规律,人们还要利用这些基本规律来研究复杂的自然现象和解决复杂的工程技术问题,这就需要求出数学物理方程中许多特定问题的解答,随着计算机的出现及计算技术的发展,即使是相当复杂的问题,也可以计算出足够精确的数值来,这对于预测自然现象的变化(如气象预报)和进行各种工程设计(如机械强度的计算)都有着很重要的作用.在研究数学物理方程的同时,人们对偏微分方程的性质也了解得越来越多、越来越深入,形成数学中的一门重要分支——偏微分方程理论.它既有悠久的历史,又不断地更新着它的对象、内容和方法.它直接联系着众多自然现象和实际问题,不断地提出或产生需要解决的新课题和新方法.它所面临的数学问题多样而复杂,不断地促进着许多相关数学分支(如泛函分析、复变函数、微分几何、计算数学等)的发展,并从它们之间引进许多有力的解决问题的工具.因此,数学物理方程又是纯粹数学的抛物型方程解的估计及其应用许多分支和自然科学各部门及工程技术等领域之间的一个重要的桥梁.2 选题背景2.1 题目类型及来源题目类型:研究论文题目来源:专题研究2.2研究目的和意义数学物理方程将数学与物理紧密地结合在一起,用精微的数学思想和方法应用于实际的物理研究中,通过物理过程建立数学模型(偏微分方程),通过求解和分析模型,对实际物理过程进一步深入理解,提出解决实际问题的途径和方法.而抛物型方程是偏微分方程中的一类,且在研究热传导、扩散等物理现象时都会遇到,具有巨大的理论价值,同时,抛物型方程的概念和性质也决定了它在工程数学,物理等方向的巨大实用价值.研究抛物型方程解的估计及其应用,有助于我们更好的理解和掌握偏微分方程理论,并在认识和了解抛物型方程广泛的使用价值的基础上,能够探索抛物型方程更广泛的应用.随着物理科学所研究的现象在广度和深度两方面的扩展,偏微分方程的应用范围更广泛.从数学自身的角度看,抛物型方程的求解促使数学在函数论、变分法、级数展开、常微分方程、代数、微分几何等各方面进行发展.2.3国内外现状和发展趋势与研究的主攻方向自18世纪以来,偏微分方程理论在得到广泛应用的同时,也得到不断发展和完善,内容也越来越丰富,它直接联系着众多自然现象和实际问题,不断地提出或产生需要解决的新课题和新方法.它所面临的数学问题多样而复杂,不断地促进着许多相关数学分支的发展,如泛函分析、复变函数、微分几何、计算数学等,并从它们之中引进许多有力的解决问题的方法.关于抛物型方程的解,有经典解、广义解、数值解等方面的研究.经典解在求解区域中具有方程中所出现的连续偏导数,并按通常意义满足方程与定解条件,也就是热传导方程的一些知识说,将解代入方程及定解条件后即可使其化为恒等式.求经典解的方法有分离变量法、Fourier变换法.经典解容易理解,且应用方便,但实际求解抛物型方程的定解问题时,往往不一定能得到经典解.于是就提出了广义解[1]的理论,即先寻求一个正则性较低的函数(广义解),它按较弱的意义满足方程与定解条件,然后再进一步证明这个函数实际上就是原来问题的经典解.广义解可按逼近过程来定义,也可按分部积分的方法来定义.由于抛物型方程的广义解是在“较差”的函数类中寻求相应定解问题按拓广意义下的解,因而又提出了广义函数的概念、性质、应用等方面的研究.在实际求解抛物型方程的定解条件时,除了一些特殊的情况下可以方便地求得其精确解外,在一般的情况下,当方程或定解条件具有比较复杂的形式,或求解区域具有比较复杂的形状时,往往求不到或不易求到其精确解,我们就只得去寻求抛物型方程定解问题的近似解,特别是数值近似解,于是就有了抛物型方程数值解的理论研究.求抛物型数值解的方法主要有:有限差分法、有限元素法.随着社会文明的发展,我国与其它国家的文化交流沟通很全面,偏微分方程的研究方向基本上也是一致的.在微分方程定性理论中有着重要应用的时滞积分不等式以及在差分方程定性理论中有重要应用的时滞离散不等式和关于时间尺度的动力方程理论,在研究了抛物型方程及抛物型方程系统、双曲型方程及双曲型方程系统的强迫振动性理论以及某些时滞脉冲偏微分方程在特定边值条件下解的振动的若干充要条件和多时滞脉冲抛物型微分方程系统解的强迫振动性,推广了已有的结果,建立了若干新定理,促进了偏泛函微分方程理论的发展,并主要创新提出并建立了偏泛函微分方程系统的强迫振动性理论以及具有连续分布变元的双曲型偏泛函微分方程系统和抛物型偏泛函微分方程系统的振动性理论;获得了变系数抛物型偏泛函微分方程解的振动的若干充分必要条件和时滞脉冲抛物型微分方程解的振动的若干充分必要条件和多时滞脉冲抛物型微分系统解的强迫振动性;研究了时滞积分不等式理论和关于时间尺度的动力不等式理论.这些理论的建立,为偏泛函微分方程理论和积分不等式理论的进一步发展起到了非常大的促进作用.3热传导方程的一些知识3.1 热传导方程的导出若物体内各点的温度不相同,则其热量就会从温度较高的地方向温度较低的地方抛物型方程解的估计及其应用传递,这就是常说的热传导现象.由于热量的传递,所以物体内温度随时间和点的位置不同而不同,因此热传导问题可归结为研究物体内部的温度分布情况.下面我们考察一个均匀的、各向同性的物体G 在Ω内部的温度变化规律. 设以(),,u x y z 表示物体G 在Ω内任一点(),,M x y z 处在时刻t 的温度.在Ω内任取一小块区域V ,使V -⊂Ω,并且其边界Γ是光滑的闭曲面,Γ上面积元素的单位外法向量记作n .根据传热学中的傅里叶实验定律[2],物体在无穷小时段dt 内,从V 内经过dS 流出的热量dQ 与时间dt ,流经面积dS 以及温度沿dS 的外法向量的方向导数un∂∂成正比,即ud Q k d S d t k u n d S d t n∂=-=-∇⋅∂ 其中0k >是物体的热传导系数,,,x y z ⎛⎫∂∂∂∇= ⎪∂∂∂⎝⎭.上式中的负号表示热流的方向与温度梯度的方向相反(因为热量总是由温度高处流向温度低处),因此从时刻1t 到时刻2t 经过Γ流入V 内的全部热量211t t Q d t k u n d s d t Γ=∇⋅⎰⎰⎰若物体Ω内有热源,且热源强度为(),,,F x y z t (即在时刻t 点(),,x y z 处的单位面积在单位时间内发出的热量),则在[]12,t t 内,V 从热源上吸收的热量为 ()212,,,t t VQ F x y z t d x d yd z d t =⎰⎰⎰⎰另一方面,在[]12,t t 内,V 内温度从()1,,,u x y z t 升高到()1,,,u x y z t 所需吸收的热量为()()321,,,,,,VQ c ux y z t u x y z td x d y d zρ=-⎡⎤⎣⎦⎰⎰⎰ 其中为c 物体的比热,ρ为物体的密度. 根据能量守恒,有热传导方程的一些知识123Q Q Q +=若(),,,u x y z t 关于,,x y z 具有二阶连续偏导数,则由高斯公式得 22111t t t t VQ dt k u ndsdt k udxdydzdt Γ=∇⋅=∆⎰⎰⎰⎰⎰⎰⎰这里 ∆ 是laplace 算子,222222x y z∂∂∂∆=++∂∂∂若(),,,u x y z t 关于t 具有一阶连续偏导数,则由Newton-Leibniz 公式有213t t t VQ d t c u d x d y d zρ=⎰⎰⎰⎰ 因此有()2211t t t t t VVdt c u dxdydz dt k u F dxdydz ρ=∇+⎰⎰⎰⎰⎰⎰⎰⎰由于时间段[]12,t t 及区域V 是任意取定的,并且被积函数是连续的,则2t u a u f-∆= 其中2k a c ρ=,Ff c ρ=,并且当0f ≥时,表示Ω内有热源;当0f ≤时,表示Ω内有冷源(即热汇).在适当情况下,方程中描述空间坐标的独立变量的数目还可以减少.例如当物体是均匀细杆时,假如它的侧面是绝热的,也就是说不产生热交换,又假定温度的分布在同一截面是相同的,则温度函数u 仅与坐标x 及时间t 有关,我们就得到一维热传导方程222u u a t x ∂∂=∂∂ 同样,如考虑薄片的热传导,薄片的侧面绝热,可得二维热传导方程22222u u u a t x y ⎛⎫∂∂∂=+ ⎪∂∂∂⎝⎭抛物型方程解的估计及其应用3.2 定解问题的提法方程描述的是同类物理现象的共性,但是每一具体的物理现象都是处在各自特定条件之下的,这就需要我们把它所处的特定条件也用数学形式表达出来,我们称这些特定条件为定解条件.定解条件分为初始条件和边界条件.初始条件是说明初始状态的条件,边界条件是描述边界状态的条件,边界条件可分为三类,第一类边界条件(又称Dirichlet 边界条件)是直接给出未知函数在研究区域Ω的边界∂Ω上的值;第二类边界条件(又称Neumann 边界条件)是在∂Ω上给出未知函数u 沿∂Ω沿外法方向n 的方向导数;第三类边界条件(又称为Robin 条件)是在边界∂Ω上给出未知函数u 及其沿∂Ω的外法方向导数的某种线性组合的值.从物理学角度来看,如果知道了物体在边界上的温度状况(或热交换状况)和物体在初始时刻的温度,就可以完全确定物体在以后时刻的温度.因此热传导方程最自然的一个定解问题就是在已给的初始条件与边界条件下求解问题的解.初始条件的提法显然为()(),,,0,,u x y z x y z ϕ=其中(),,x y z ϕ为已知函数,表示物体在0t =时的温度分布第一边界条件:在3R 中的有界区域Ω的导热问题中,若Ω的边界∂Ω处于恒温0u 的环境下,则边界条件为0u u ∂Ω|=若边界温度按已知规律(),,,g x y z t 变化,则(),,,u g xy z t ∂Ω|= 第二边界条件:若热量在边界曲面∂Ω各点的流速为(),,,G x y z t ,则由Fourier 定律,边界条件可写成(),,,ug x y z t n ∂=∂ 其中Gg k =-,若0G =,则0u n ∂Ω∂=∂,此时称之为绝热边界条件.定解问题的求解第三边界条件:如果物体内部与周围的介质通过边界∂Ω有热量交换,物体外介质的温度为2u ,物体表面的温度为1u ,内外两种介质间的热交换系数为()110k k >,根据Newton 定律,从物体内部流到外部的热量与两介质间的温度差成正比,即有()112d Q k u u d sd t =-另一方面,由Fourier 定律[3],在时间间隔内从边界曲面上面积元流出的热量为ud Q k d s d tn∂=-∂ 从而有()112u k u u d s dt k d s d t n∂-=-∂ 即(),,,u u g x y z t n σ∂Ω∂⎛⎫+=⎪∂⎝⎭ 其中1k kσ=, ()1,,,u g x y z t σ= 4 定解问题的求解4.1 初值问题的求解我们可以利用傅里叶变换来求解热传导问题的初值问题.其思想是把原函数变换到另一类函数中去,经过变换,使热传导方程变为常微分方程,从而可以找出一个解,再经过Fourier 的逆变换,得到原热传导方程的解.()()()()2,,,,0,t xx yy u a u u f x y t u x y x y ϕ⎧-+=⎪⎨=⎪⎩ (1) 视t 为参数,先求解齐次热传导方程的初值问题()()()2,,0,t x x yy u a u u u xy x y ϕ⎧=+⎪⎨=⎪⎩ (2)对,x y 进行Fourier 变换,记()()12,,,,F u x y t U t λλ=⎡⎤⎣⎦,抛物型方程解的估计及其应用()()12,,F x y ϕλλ=Φ⎡⎤⎣⎦在(1)式两边关于,x y 进行Fourier 变换,原问题变为()()()()()()()222121122121212,,,,,,,,0,d U t a i U t i U t dtU λλλλλλλλλλλλ⎧⎡⎤=+⎪⎣⎦⎨⎪=Φ⎩(3) (2)式是带参数12,λλ的常微分方程的柯西问题,它的解为()()()2121212,,,a tU t eλλλλλλ-+=Φ (4)函数()212a teλλ-+的Fourier 逆变换[4]为()()()()()()()2222221212122222112211221221F 21=2a a i x y a t i xa t i xe t e te d d ed ed λλλλλλλλλλλλπλλπ+∞-+-++-∞----+∞+∞-∞-∞⎡⎤=⎢⎥⎣⎦⎰⎰⎰-()222222111122111111+11cos sin =2cos a t i xa ta ta ted exd i exd exd λλλλλλλλλλλλ----+∞+∞+∞-∞-∞-∞∞-=+⎰⎰⎰⎰令()221+110cos a t I x e xd λλλ∞-=⎰()()221222211+/111111202sin 1 =sin cos 2 =2a t a t a t Ix e xd e x x xe d a t xI x a tλλλλλλλλλ∞-+∞--+∞0=-⎡⎤∣-⎢⎥⎣⎦-⎰⎰ 解得()224x a tI x ce-=又()2212+1+00 a t y I e d e dy λλ∞-∞-===⎰⎰定解问题的求解则有()222222121421F 4x y a a tet e a tλλπ+--+⎡⎤=⎢⎥⎣⎦-由(4)可得初值问题(2)的解为()()()()222421,,,4x y a tu x y t ed d a t ξηϕξηξηπ-+--+∞+∞-∞-∞=⎰⎰ (5)再求解非齐次热传导方程具有齐次初始条件的柯西问题()()()2,,,.00t xx yy u a u u f x y t u x y ⎧=++⎪⎨=⎪⎩ (6) 由齐次化原理[5],此柯西问题的解可写为()(),,,,;tu xy t x y t d ωττ=⎰ 而(),,;x y t ωωτ=为下述柯西问题的解:()()()2,,,,,t x x yy a t x y f x y ωωωτωττ⎧=+>⎪⎨=⎪⎩于是,利用(5)式,易知柯西问题(6)的解为()()()()()222420,,1,,4x y t a t u x y t ed d d a t ξητϕξητξητπτ-+--+∞+∞--∞-∞=-⎰⎰⎰ (7)由叠加原理[6],由(5)及(7)就得到柯西问题(1)的解为()()()()()()()()222222424201,,,4,,1 4x y a tx y ta t u x y t ed d a t ed d d a t ξηξητϕξηξηπϕξητξητπτ-+--+∞+∞-∞-∞-+--+∞+∞--∞-∞=+-⎰⎰⎰⎰⎰在上面的推导中,由于预先不知道是否满足进行傅里叶变换及有关计算的条件,所得的解还只是形式解.为证明上式确实是柯西问题(1)的解,还得进行验证.抛物型方程解的估计及其应用4.2 初边值问题的求解热传导方程的初边值问题20 t xx u a u -= (8)00x x l u u ==∣=∣= (9) ()0 t u x ϕ=∣= (10) 令()()() ,u x t X x T t = (11)并要求它满足齐次边界条件(9),这里()X x 及()T t 分别表示仅与x 有关及仅与t 有关的特定函数.将(11)代入方程(8)中,得到()()()()///X x T t X x T t -= (12) 将上式分离变量,有()()()()///2T t X x a T t X x λ==- (13) 由于在(13)式中,左边仅是t 的函数,右边仅是x 的函数,左右两端要相等,只有等于同一个常数才可能.记次常数为λ-(其值待定),就得到()()/2T aT 0t t λ+= (14)()()//0Xx X x λ+= (15)这样方程(13)就被分离为两个常微分方程,其中一个含有自变量t ,另一个仅含有自变量x ,我们可以通过求解这两个方程来决定()T t 及()X x ,从而得到方程(8)的特解(11)为了使此解是满足齐次边界条件(9)的非平凡解,就必须找到方程(8)满足边界条件定解问题的求解()()00,0X X l == (16) 的非平凡解.方程(15)的通解随0λ>,0λ=以及0λ<而不同,下面分三种情况讨论:情形1 当0λ<时,方程(15)的通解可写成 ()12X x C C e =+要使它满足边界条件(16),就必须1200C C e +=⎧⎪⎨+=⎪⎩ 由于110e≠只能120C C ==.故在0λ<的情况得不到非平凡解.情形2 当0λ=时,方程(15)的通解可以写成 ()12X x C C X =+ 要满足边界条件(16),()X x 也只能恒等于零.情形3 当0λ>时,方程(15)的通解具有如下形式: ()12X x C C =+ 由边界条件()00X =知10C =,再由()2s i 0X l C l ==可知,为了使20C ≠,就必须0=.于是222(1,2,)k k k lπλλ===⋯这样就找到了一族非零解()s i n(1,2,)k k k X x A x k lπ==⋯ 将固有值代入方程(14)中,可得到其通解为()2222(1,2,)a k tl k k T t B ek π-==⋯ 这样就得到方程(8)的满足齐次边界(9)的下列分离变量形式的特解:()()()2222k ,s i n (1,2,)a k tl k k k k u x t X x T t a ex k lππ-===⋯现在我们设法作这种特解的适当的线性组合,以得出初边值问题的解,也就是说,要决定常数k a 使()22221,s i n a k tl k k k u x t a ex lππ∞-==∑ (17)满足初始条件(10). 故由初始条件(10)应有()1s i n kk k x a x l πϕ∞==∑ 由于 1,sin k x l π⎧⎫⎨⎬⎩⎭在[]0,l 上正交,因此,k a 是在[]0,l 区间中正弦展开的傅里叶级数的系数,即()02sin l k k a d l lπϕξξξ=⎰ (18) 故()()222201,sin sina k tll k k k u x t d ex l lπππϕξξξ∞-==⋅∑⎰ (19) 是用级数形式表示的初边值问题的形式解.为了考察由分离变量法得到的形式解是否是混合问题的经典解,还得进行验证. 当1C ϕ∈,且()()00l ϕϕ==,()x ϕ是有界函数,(18)式确定的函数(),u x t 是混合问题的解.分析:在求解过程中,级数(17)中的每一项都满足方程(8),因此只要证明级数(17)可以逐项求导两次就好了.也就是说,如果证明了级数(17)求导两次后仍是一致收敛的,那么它一定满足方程(8),此时边界条件(9)和初始条件(10)的满足也是显然的推论了.证明:由于式(19)中含有因子2222a k tl eπ-,因此对于任意0δ>,当0t >时,对任意的0p >,级数22221p a k tl k k el ππ∞-=⎛⎫⎪⎝⎭∑均是一致收敛的,而由ϕ是有界函数的假设(()x M ϕ<),可得()0sinlk d Ml lπϕξξξ≤⎰故(19)式中列举的所有级数是一致收敛的,因而,由式(19)表示的级数,当0t >时,关于x 及t 是无穷次可导的,并且求导与求和可以交换.由于级数的每一项都满足方程(8)及边界条件(9)、(10),从而式(19)式表示的级数在0t >时确实满足方程及边界条件.当加上条件()()00l ϕϕ==时,当0t →时,对任意[]0,x l ∈,由式(19)给出的级数趋于初值()x ϕ,即得到式(19)给出的级数确实是初边值问题(8)~(10)的经典解.5 抛物型方程解的估计及其应用先验估计是偏微分方程理论研究中的一个常用的方法.其特点是在假设定解问题解存在的前提下导出解所应当满足的估计,而常用的估计有最大模估计[7],能量估计[8]等等.一般地,我们可以根据先验估计得到定解问题解的唯一性和稳定性,并且可结合其他一些分析方法推导出解的存在性,此外,作为对解的一种估计,先验估计还可能提供关于解的某种性态(如有界性等)方面的信息.5.1 极值原理考虑热传导方程()()2,,,t x x L u u a u f x t x t Q≡-=∈ 其中(){},0,0Q x t x l t T =<<<≤,Q 的侧边和底边统称为Q 的抛物边界,记作Γ,即(){}(){}(){},0,0,,0,0,0x t x t T x t x l t T x t t x l Γ==<≤⋃=<≤⋃=≤≤在热传导过程中,如果物体内部无热源,则热量总是由温度高处向其它地方扩散,而温度最低处的温度会逐渐上升.因此物体的最高温度和最低温度总是在初始时刻或物体的边界上达到.这就是热传导方程的“极值原理”.定理 1(弱极值原理) 设函数()()()2,1,C u x t Q C Q ∈⋂满足Lu f =. (1) 若0f ≤,则u 在Q 上的最大值必在抛物边界Γ上达到,即 ()()m a x ,m a x ,Qu x t u x t Γ=(2) 若0f ≥,则()()m i n ,m i n ,Qu x t u x t Γ=(3) 若0f =,则()()m a x ,m a x ,Qu x t u x t Γ=, ()()m i n ,m i n ,Qu x t u x t Γ=同时成立,这里()2,1C Q 表示在Q 内关于x 二次连续可微,且关于t 一次连续可微的函数全体.证明:(1)不妨先考虑0f <情形. 反设存在点()00,x t Q ∈,使得()()00,max ,Qu x t u x t =则在该点处0x u =,0xx u ≤,0t u ≥(如果0t T <,则0t u =;如果0t T =,则0t u ≥).因此()()()00200,,0t xx x t f x t u a u =-≥,这与0f <的假设相矛盾.故(),u x t 不能在Q 内达到最大值,从而有()()m a x,m a x ,Qu x t u x t Γ= 当 (),0f x t ≤时,设法将其转化为前面的情形.为此构造辅助函数 ()(),,v x t u x t t ε=- 其中ε是任意小的正数.因为0L v L u f εε=-=-<所以()()m a x ,m a x ,Qv x t v x t Γ=于是()()()()max ,max ,max ,max ,QQu x t v x t t v x t T u x t T εεεΓΓ=+≤+≤+⎡⎤⎣⎦令0ε→,得()()m a x ,m a x ,Qu x t u x t Γ=(2)若0f ≥,则对u -应用情形(1)的结论即可.(3)结合前面两种情况,若0Lu =,则u 在Q 的上的最大值与最小值都在抛物边界Γ上达到.下面我们将弱极值原理推广到稍一般的热传导方程()()()21,,,t x x x L u u a u bx t u c x t u f x t≡-++= 定理 2 函数()()2,1u C Q C Q ∈⋂满足10L u f =≤,则u 在Q 上的正最大值必在抛物边界Γ上达到,即()()m a x ,m a x ,Qu x t u x t +Γ≤由于其证明与定理1的证明方式类似,这里不再赘述.定理3 设()0,c x t c ≥-,其中0c 为正常数.若函数()()()2,1,u x t C Q C Q ∈⋂满足10L u f =≤,且()max ,0u x t Γ≤,则必有()max ,0Qu x t ≤证明 令()()0,,c t v x t e u x t -=,则(),v x t 满足方程 ()0200c t t xx x v a v bv c c v fe --+++=≤ 由于00c c +≥,根据定理2,得()()()0m a x,m a x ,m a x ,0c tQv x t v x t e u x t -++ΓΓ≤≤≤ 因此结论得证.利用定理3,不难得到下列推论:推论1(比较原理) 设()()00,0c x t c c ≥-≥,又设()()2,1,u v C Q C Q ∈⋂,且11L u L v ≤,u v ΓΓ≤,则对任意的(),x t Q ∈,有 ()(),,u x t v x t ≤5.2 初边值问题解的最大模估计设Ω是n R 中的有界开集,0T >.记(0,]T Q T =Ω⨯,(){}()[0,)0T T Γ=∂Ω⨯⋃Ω⨯这里的T Γ称为T Q 的抛物边界.我们先在T Q 中研究抛物型方程记 []()()1,,int i x i A u u u b x t uf x t ==-∆+=∑[]()()()1,,,i nt ix i B u u u b xt u c x t u f x t==-∆++=∑ 考察第一初边值问题[]()()()()()()()()[]1,, ,,0 ,, ,0,i nt i x Ti A u u u b x t u f x t x t Qu x x xu x t g x t x t T ϕ=⎧=-∆+=∈⎪⎪⎪=∈Ω⎨⎪=∈∂Ω⨯⎪⎪⎩∑ (20)定理4 设()()2,1T T u C Q C Q ∈⋂是问题(20)的解,则TQ max u FT B ≤+其中sup TQ F f =,()[]{}0,max max ,max T B x g ϕ∂Ω⨯Ω=证明 令v tF B =+,与u ±作比较.因为 [][]A u F f A u =≥±=± ,(),T x t Q ∈()()(),0,0v x B x u x ϕ=≥±=± , x ∈Ω v B g u ∂Ω∂Ω∂Ω≥≥±=± , 0t T ≤≤ 由比较原理知,u v FT B ±≤≤+,即 ()TQ max ,u x t FT B ≤+推论 2 第一初边值问题(20)的解在函数类()()2,1T T C Q C Q ⋂中是唯一的,且连续地依赖于f ,ϕ和g .证明 当0f g ϕ==≡时,对应的解u 满足TQ max 0u =,故0u ≡,从而解是唯一的.假设i u 是对应于{},,i i i f g ϕ的解,1,2i =,则12u u -是对应于{}121212,,f f g g ϕϕ---的解.于是[]{}TT121212120,Q Q max max ax max ,max T u u T f f g g ϕϕ∂Ω⨯Ω-≤-+--所以当{}111,,f g ϕ与{}222,,f g ϕ充分接近时,1u 与2u 也充分接近,这说明问题(20)的解连续地依赖于f ,ϕ和g .现在考察第一初边值问题[]()()()()()()()[], ,,0 ,, ,0,TB u f x t x t Q u x x x u x t g x t x t T ϕ⎧=∈⎪⎪=∈Ω⎨⎪=∈∂Ω⨯⎪⎩ (21) 定理5 设()0,c x t c ≥-,()()2,1T T u C Q C Q ∈⋂是问题(21)的解,则 ()0TQ m a x c T u e FT B ≤+其中sup TQ F f =,()[]{}0,max max ,max T B x g ϕ∂Ω⨯Ω=.证明 不妨认为00c ≥,令()0c t v e FT B =+,与u ±作比较.因为[]()()()()()()[]()00000000, =, ,c t c t c t c t c t c t T B u Fe c e Ft B c x t e Ft B Fe e c c x t Ft B Fe F f B u x t Q =+++++++≥≥≥±=±∈()()(),0,0v x B x u x ϕ=≥±=± , x ∈Ω v B g u ∂Ω∂Ω∂Ω≥≥±=± , 0t T ≤≤ 由比较原理知,()0c Tu v e FT B ±≤≤+,即()()0TQmax, c T u x t e FT B ≤+5.3 初值问题解的最大模估计记[]T D 0,n R T =⨯,[](),t C u u u c x t u =-∆+ 考察初值问题[]()()()(), ,,0 TnC u f x t x tD u x x x Rϕ⎧=∈⎪⎨=∈⎪⎩ (22) 设(),c x t 连续,()()00,0c x t c c ≥->,(),f x t 和()x ϕ有界,记 s u p TD F f =, sup nR ϕΦ=如果()()2,1T T u C D C D ∈⋂是初值问题(22)的解,则 ()0s u p Tc T D u e FT ≤+Φ证明 令()()0,,c t v x t u x t e -=,则v 满足[]()()()(),,,0 t nD v v v c x t v f x t v x x x R ϕ⎧=-∆+=⎪⎨=∈⎪⎩ (23) 其中()()0,,0c x t c x t c =+≥,()()0,,c t f x t e f x t -=由于解得先验估计方法不能直接用于初值问题,我们希望借助于一个有界区域上的初边值问题进行讨论,任意取定较大的常数L ,记{}](,0,L T D x L T =≤⨯.因为解u 有界,所以存在正常数K 使得u K ≤在D T 上成立,在有界区域,L T D 上考虑辅助函数()()22,2K w x t Ft x nt v L =+Φ++± 直接计算知,在,L T D 上w 满足[]()()()()()()002,22220 ,,0 ,,0c t L T c tx L x L K D w F c Ft x nt e f x t D L K w x x x x LL w x t K u x t e ϕ--==⎧⎧⎫=++Φ++±≥∈⎨⎬⎪⎩⎭⎪⎪=Φ+±≤⎨⎪⎪≥Φ+±>⎪⎩利用比较原理知,(),0w x t ≥在,L T D 上成立对于D T 内的任一点()00,x t ,取L 充分大使得()00,,L T x t D ∈,于是()00,0w x t ≥ 即()()2000002,2K v x t Ft x nt L≤+Φ++ 令L →∞得()000,v x t Ft Ft ≤+Φ≤+Φ从而()()()000000,,c t c T u x t v x t e e Ft =≤+Φ由()00,T x t D ∈的任意性知,估计式(23)成立.推论3 初值问题(23)的解在函数类()()2,1T T C D C D ⋂中是唯一的,且连续地依赖于f ,ϕ.由于其证明与推论3的证明方式类似,这里不再赘述.5.4 初边值问题的能量估计设Ω是n R 中的一个光滑区域,在](0,T Q T =Ω⨯上考察第一初边值问题()()()()()[], ,,0 0 ,0,t T u u f x t x t Q u x x x u x t T ϕ-∆=∈⎧⎪⎪=∈Ω⎨⎪=∈∂Ω⨯⎪⎩ (24) 定理6 设()()1,02,1T T u C Q C Q ∈⋂是问题(23)的解,则存在正常数()C C T =使得()222200max ,2TT t Tux t dx u dxdt C dx f dxdt ϕΩΩΩΩ≤≤⎛⎫+∇≤+ ⎪⎝⎭⎰⎰⎰⎰⎰⎰ (25) 证明 问题(24的方程两边乘以u 并在T Q 上积分,得000tttt uu dxdt u udxdt f udxdt ΩΩΩ-∆=⎰⎰⎰⎰⎰⎰(26)对(26)式左端第一项中关于t 的积分利用分部积分以及初值条件,可知()()22011,22t t uu dt u x t x ϕ=-⎰ (27)对(26)式左端第二项关于x 的积分利用散度定理以及边界条件,推出22u u u d x u d S u d x u d x n Ω∂ΩΩΩ∂∆=-∇=-∇∂⎰⎰⎰⎰ (28) 将(27)式和(28)式代入(26)式,得2220022ttu dx u dxdt f udxdt dx ϕΩΩΩΩ+∇=+⎰⎰⎰⎰⎰⎰ (29)利用不等式222ab a b ≤+可知 220002t ttf u d x d t f d x d tud x d t ΩΩΩ≤+⎰⎰⎰⎰⎰⎰ 将上式代入(29)式,得222220002tttu dx u dxdt f dxd u dxdt dx ϕΩΩΩΩΩ+∇≤++⎰⎰⎰⎰⎰⎰⎰⎰ (30)记 ()20tY t u dxdt Ω=⎰⎰,()220t F t f dxd dx ϕΩΩ=+⎰⎰⎰那么不等式蕴含()()()Y t Y t F t '≤+ 利用Gronwall 不等式[9]推出()()()()()()2022001 tt t t ttu dxdtF t Y t Y e e F t e F t e f dxd dx ϕΩΩΩ=≤+-⎛⎫≤=+ ⎪⎝⎭⎰⎰⎰⎰⎰将上式代入(30)式知()22220021tt tu dx u dxdt e f dxd dx ϕΩΩΩΩ⎛⎫+∇≤++ ⎪⎝⎭⎰⎰⎰⎰⎰⎰ 此式两边关于t 在[]0,T 上取上确界,就得到估计式(25).下面我们将讨论一般形式的二阶抛物型初边值问题.设Ω为n R 中的有界区域,且有光滑边界()0,T Q T =Ω⨯,在区域中讨论一般形式的二阶抛物型初边值问题()()()(),11,,,,i j i nnij x x i x i j i u a x t u u b x t u c x t u f x t t ==∂-++=∂∑∑ (31) ()0 t u x x ϕ==∈Ω (32) 0T u ∑= (33) 解的性质.式中,()0,T T ∑=Γ⨯为区域的侧边界;()12,,n x x x x =∈Ω为方便讨论,作如下假设:(1) 系数ij a 、i b 、c 及右端项f 都是T Q 上的连续函数,并且ij a 在T Q 上还具有一阶连续偏导数.(2) 对一切,1,2,i j n = ;ij ji a a =且存在正常数0α>,使得对一切(),T x t Q ∈及任意给定的实向量()12,,,n ξξξ ,有:()2,11,nnijiji i j i a x t ξξαξ==≥∑∑成立.对于初边值问题的解,定义能量函数:()212E t u d x Ω=⎰ (34)定理7 若(),u x t 为初边值问题(31)~(33)的解,能量函数()E t 按式(34)定义,则能量估计式:()()200 0t C tC tE t E e C e f d x d t t T Ω≤+≤≤⎰⎰(35)成立.其中,C 为一个不依赖于u 的正常数.证明 用u 乘以式(31),并在Ω上关于x 积分,就得到:()(),11,,i j i n n t ij x x i x i j i u udx a x t u u dx b x t u u cu dx fudx ΩΩΩΩ==⎛⎫⎛⎫-++= ⎪ ⎪⎝⎭⎝⎭∑∑⎰⎰⎰⎰ []0,t T ∈ (36) 式左端的第一项可以写成212d u dx dt Ω⎧⎫⎨⎬⎩⎭⎰;当3n ≥时,记12,,,n ααα 为侧边界T ∑法向量的方向角,dS 为广义面积微元.令(),1,2,,i ij ij x p a uu i j n == ,固定i ,让1,2,,j n = ,利用高维高斯公式[10],并注意边界条件(它隐含着0T u ∑=),边界积分项为零,可得()()()()()()12121122121212120cos cos cos = ==Ti i i n i ni i in n i i in n i x x i x x in x x i i in x x x x p p p dSp p p dx x x x a u a u a u udx a u a u a u u dxααα∑ΩΩΩ=++⎛⎫∂∂∂+++ ⎪∂∂∂⎝⎭+++++⎰⎰⎰⎰故对固定的i ,有:()()()()()12121212=i i i n i ni x x i x x in x x i i in x x x x a u a u a u udx a u a u a u u dx ΩΩ-+++++⎰⎰(37)成立,对式(37)关于i 从1到n 求和.式(36)左端的第二项可以写成:(),1,1,1i j i j iin n ni j x x i j x x i jx x i j i j i j a u u d x a u u d x au u d x ΩΩΩ===⎛⎫⎛⎫⎛⎫-=+ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭∑∑∑⎰⎰⎰ (38) 将上式的第二项,连同式右端的第三、四项移至等式右边,并将其和记为(),t Q u u dx Ω⎰ 则有()()()1,1,,i i i n n ti x ij x x i i j Q u u dx fudx b x t u u cu a u u dx ΩΩΩ==⎛⎫=-++ ⎪⎝⎭∑∑⎰⎰⎰则由于系数的可微性假设(1)可得,对一切0t T ≤≤成立()21,i nt T xi Q u u d x C u u u d x ΩΩ=⎛⎫≤+ ⎪⎝⎭∑⎰⎰ (39) 其中T C 为一个不依赖于T 的正常数,但与u 无关.对任意给定的0ε>,有2211122i innnx x i i i nuu dx u dx u dx εεΩΩΩ===≤+∑∑∑⎰⎰⎰ (40)取TC αε=,由式(40)就得到()22111,2innt x i i Q u u dx u dx C u dx αΩΩΩ==≤+∑∑⎰⎰⎰(41)其中212TT nC C C α=+,将式(41)代入式(36),容易得到2221,11111222i j i n n n ij x x x i j i i dE a u u dx u dx C u dx f dx dt αΩΩΩΩ===⎛⎫⎛⎫+≤+++ ⎪ ⎪⎝⎭⎝⎭∑∑∑⎰⎰⎰⎰ (42) 再注意到由假设(2)有2,11i j i n n ij x x x i j i a u u dx u dx αΩΩ==⎛⎫≥ ⎪⎝⎭∑∑⎰⎰ 就可得到()22dEC E t f dx dtΩ≤+⎰ (43)其中2121C C =+.在式(43)两边乘以2C t e -再对t 积分,,并放大被积函数,即可得 ()()200tC tC tE t E e C ef dx d t Ω≤+⎰⎰定理证毕.5.5 能量不等式的应用5.5.1 初边值问题解的唯一性热传导方程是抛物型方程的典型代表.下面考虑二维热传导方程的初边值问题()2t xx yy u a u u f =++ (44)()0,t u x y ϕ== (45) (),,u x y t μΓ= (46) 这里,Γ表示Ω的边界,应用能量不等式可得如下定理.定理8 若热传导方程的初边值问题的解存在,则其解唯一.证明 设1u ,2u 是该定解问题的两个解,则其差12u u u =-满足相应的齐次方程及齐次初始条件和齐次边界条件.此时的齐次方程满足假设(1)、(2),有(34)式定义的能量函数知,在初始时刻有()00E =,故由能量不等式(35)得:()()22220x y E t u a u u dxdy Ω⎡⎤=++=⎣⎦⎰⎰ 即0x y u u u ===,从而可推出(),,u x y t const =.又由于在初始时刻0u =,故得(),,0u x y t ≡.即12u u =.这样就证明了初边值问题(44)~(46)解的唯一性. 5.5.2 初边值问题解的稳定性为了记号简单起见,对于定义在区域Ω上的函数ϕ和定义在区域上()0,T ⨯Ω的函数f ,常以()2L ϕΩ和()()20,L T f ⨯Ω分别表示()122dxdy ϕΩ⎰⎰和()1220Tf dxdydt Ω⎰⎰⎰.定理9 热传导方程的初边值问题:()2t xx yy u a u u f =++()0,t u x y ϕ== 0u Γ=的解(),,u x y t ,在下述意义下关于初始值ϕ与方程右端项f 是稳定的:对任何给定的0ε>,一定可以找到仅依赖于ε和T 的0η>,只要 ()212L ϕϕηΩ-≤ ()212x xL ϕϕηΩ-≤()212y yL ϕϕηΩ-≤ ()()2120,L T f f η⨯Ω-≤ (47)那么以1ϕ为初值、1f 为右端项的解1u 与以2ϕ为初值、2f 为右端项的解2u 之差在上满足()212L u u εΩ-≤ ()212x xL u u εΩ-≤ ()212y yL u u εΩ-≤ (48)证明 记12u u u =-,12ϕϕϕ=-,1f f f =-,则u 满足()2t xx yy u a u u f =++ (49)()0,t u x y ϕ== (50) 0u Γ= (51) 方程(49)满足假设(1)、(2),从而利用能量不等式(35),可得:()()()()222000t TCt Ct E t E e Ce f dxdydt C E f dxdydt ΩΩ≤+≤+⎰⎰⎰⎰⎰⎰[]0,t T ∈ (52)式中,2C 为一个仅依赖于T 的正常数.记。
§6非定常流的有限元方法考虑一个 2m 阶的抛物型微分方程定解问题(初边值问题): 求 ()()2,m u C 孜W x ,()()1,0,u t C T 轾孜犏臌,满足 112220 , 00 , 0, 0 t u A u f t t B u t B u g t u f =ìï¶ï+=W >ïï¶ïïïïï=G ïïíïïï=G ïïïïïï=Wïïïî(6.1)式中 A 仍是一个 2m 阶椭圆型空间微分算子。
根据 §3中的分析,只需考虑齐次强制边界条件。
我们还是像 §3那样,将它改写成积分形式(变分方程)。
首先,任取函数 ()0v C ¥蜽 ,用它乘以方程的两边,并在 W 上积分,得uvd Auvd fvd t WWW¶W+W=W ¶蝌上式除了左边第一项之外,其余两项仍可以按照 §3中的步骤改写成某个Sobolev 空间 V 上的变分方程()(),B u v F v = , v V "原有的第一项,可将时间导数写到空间积分外面,即uvd uvd t t WW骣抖÷ç÷çW=W ÷ç÷抖ç桫蝌并推广到可积函数空间(广义的,指Lebesgue 可积)()20,L T 轾犏臌 上。
还有,定解条件中的初始条件也要写成积分形式()0t u vd vd f =WWW=W 蝌这些加在一起,就得到非定常问题的变分方程:求 u U Î ,满足()()()()()0,,,t u v B u v F v t u v v F =ìï¶ï+=ïï¶ïíïïï=ïïî, v V " (6.2)这里,空间()()()(){}2, , , ,0, U u x t u V u t L T 轾=孜孜犏臌x 内积 (),u v 表示积分uvd WW ò ,(),B u v 和 ()F v 仍是空间 V上的双线性泛函和线性泛函,初值泛函 ()v vd F f W=W ò。
一类二阶二维常系数抛物型方程的参数估计方法抛物型方程是一类很常见的偏微分方程,在数学和工程领域都有着广泛的应用。
通常情况下,我们可以通过已知的数据和条件来求解这类方程的解析解,但是有时候我们并不能得到方程的解析解,这时我们就需要考虑用参数估计的方法来求解方程的解。
在本文中,我们将讨论一类二阶二维常系数抛物型方程的参数估计方法。
这类方程的一般形式可以写为:$$\frac{∂^{2} u}{∂t^{2}} = a\frac{∂^{2} u}{∂x^{2}} +b\frac{∂u}{∂x} + cu$$其中$a,b,c$是常数,$u(x,t)$是变量。
接下来我们将介绍一种基于最小二乘法的参数估计方法,该方法可以帮助我们估计方程中的未知参数。
首先,我们假设我们有一些已知的数据点$(x_{i},t_{i},u_{i}),i=1,2,...,n$,我们可以用数据点构建一个误差函数,使其最小化。
误差函数的一种形式可以写为:$$E(a,b,c) = \sum_{i=1}^{n} (u_{i} - u(x_{i},t_{i}))^{2}$$其中$u(x_{i},t_{i})$是抛物型方程的解析解在点$(x_{i},t_{i})$的取值。
我们的目标是找到$a,b,c$的最优估计值,使得误差函数最小化。
接下来,我们可以用最小二乘法来估计$a,b,c$的值。
最小二乘法的基本思想是通过最小化误差函数来估计参数值。
我们可以将误差函数$E(a,b,c)$对$a,b,c$分别求偏导,并让导数为0,得到方程组:$$\begin{cases}\frac{∂E}{∂a} = 0 \\\frac{∂E}{∂b} = 0 \\\frac{∂E}{∂c} = 0 \\\end{cases}$$解这个方程组,我们就可以得到$a,b,c$的最优估计值。
最后,我们可以用已知的数据点和估计值$a,b,c$来重新求解抛物型方程,得到近似解,并且可以用残差分析方法来检验我们的估计结果是否准确可靠。
有相变的变系数抛物型微分方程自由边界问题有相变的变系数抛物型微分方程自由边界问题是数学分析中一个具有重要理论意义和实际应用价值的问题。
以一维空间有相变系数抛物型微分方程为例,本文分析了论文中有关自由边界问题的解析解,研究了系数变化带来的困境,具体来说,讨论了系数跳跃和变动范围的情况,提出了一称为补偿技术的解决方案,并给出了实际应用的实例,为此问题提供了有效的数值计算方法。
一维有相变系数抛物型微分方程一般形式为:$$begin{cases}frac{d^2 u}{dt^2} + q(t) frac{d u}{dt} + p(t) u(t) = f(t) t in [t_0,t_1]u(t_0) = alphau(t_1)= betaend{cases}$$其中$u(t)$为未知函数,$p(t)$、$q(t)$、$f(t)$为带有变化情况的给定函数,$t_0$、$t_1$为区间的端点, $alpha$、$beta$为初始和终止条件。
考虑系数 $p(t)$ $q(t)$在$[t_0,t_1]$内发生跳跃或变动,出现系数变化的情况。
首先,在此情况下,上述方程组可以转化为一组以多个子区间分割的子方程组:$$begin{cases}frac{d^2 u}{dt^2} + q_i(t) frac{d u}{dt} + p_i(t) u(t) = f_i(t) t in [t_{i-1},t_i]u(t_{i-1}) = u_{i-1}u(t_i)= u_iend{cases}$$其中 $i =1,2,…,n$,$n$为区间的长度,$u_{i-1}$和$u_i$分别是子方程在$[t_{i-1},t_i]$间上的初值和终值,$p_i(t)$ 、$q_i(t)$和$f_i(t)$分别代表跳跃之后的系数。
由于系数发生变动,困难就出现了:子方程组的右端项$f_i(t)$和边界条件$u_{i-1}$和$u_i$可能与系数变动之前有所差别,但是在系数变动之后,需要将这些右端项和边界条件补偿,以保证其可以与系数变动之后的子方程组相协调。
数学中的抛物型方程抛物型方程(parabolic equation)是数学中一类重要的偏微分方程,它在物理学、工程学和社会科学等领域中具有广泛的应用。
本文将从抛物型方程的定义、特征和解法等方面进行论述,以帮助读者更好地理解和应用抛物型方程。
一、抛物型方程的定义在数学中,抛物型方程是一类二维或三维偏微分方程,其形式可以表示为:∂u/∂t = a∇²u + bu + c其中,∂u/∂t 表示函数 u 对时间 t 的偏导数,∇²u 表示函数 u 对空间坐标的拉普拉斯算子,a、b、c 是常数。
抛物型方程通常描述了某一物理现象随时间变化的规律,比如热传导、扩散等。
通过解抛物型方程,我们可以预测和分析这些物理现象。
二、抛物型方程的特征1. 热传导方程抛物型方程在热传导方程中的应用是最常见的。
热传导方程描述了物体内部温度随时间和空间的变化情况。
在一维情况下,热传导方程具有以下形式:∂u/∂t = α∂²u/∂x²其中,u(x, t) 表示在时刻 t 位置为 x 的温度,α 是热扩散系数。
2. 扩散方程抛物型方程在扩散方程中的应用也是非常重要的。
扩散方程描述了物质在浓度梯度驱动下的扩散过程。
在一维情况下,扩散方程具有以下形式:∂u/∂t = D∂²u/∂x²其中,u(x, t) 表示在时刻 t 位置为 x 的物质浓度,D 是扩散系数。
三、抛物型方程的解法对于抛物型方程,我们通常采用偏微分方程的求解方法,如分离变量法、格林函数法等。
1. 分离变量法分离变量法是一种常用的求解抛物型方程的方法。
它的基本思想是将多元函数分解为几个一元函数的乘积,并利用分离后的一元函数满足各自的方程来求解。
以热传导方程为例,我们可以将其分离变量为时间部分和空间部分:u(x, t) = X(x)T(t)代入原方程,得到两个方程:X''(x)T(t)/X(x) = T'(t)/T(t) = -λ²其中,λ² 是常数。
二维变系数抛物型方程的一个高阶ADI差分格式马小霞;颜晓琳;陈汝栋【摘要】A high accuracy alternation direction implicit scheme (ADI) for solving the two-dimensional parabolic equations is presented, and the scheme is absolutely stable and the truncation error is O(τ2+h4). The experiments show the scheme is effective and advantage, and the theory is right by a numerical example.%针对二维变系数抛物型方程,构造出了一个高精度、恒稳定的交替方向隐式(ADI)差分格式,格式的截断误差阶达O(τ2+h4)。
通过数值实验,验证了理论分析的正确性和差分格式的精确性与有效性。
【期刊名称】《天津工业大学学报》【年(卷),期】2014(000)001【总页数】5页(P77-80,84)【关键词】抛物型方程;ADI格式;截断误差;恒稳定【作者】马小霞;颜晓琳;陈汝栋【作者单位】焦作大学基础部,河南焦作 454003;天津工业大学理学院,天津300387;天津工业大学理学院,天津 300387【正文语种】中文【中图分类】O241.82抛物型方程在处理废料污染、渗透、驱动、海水入侵以及半导体等工程实际问题中有着广泛的应用,因此研究其高精度、高稳定和计算量较小的数值解法具有重要的意义.用有限差分方法研究这类问题的数值方法目前已做了许多工作[1-5].但这些工作大多是对常系数而言的.文献[4]中对二维变系数抛物型方程数值方法仅对系数依赖于一个变量的情况进行了研究,本文的研究是对系数依赖于两个变量的情形进行的.应用Taylor展开、算子方法[6]以及粘结系数法[7]得到了一个高精度(截断误差阶达O(τ2+h4))、恒稳定的ADI格式.格式的建立和稳定性分析都比文献[4]简单得多,文末的数值实验证明了本文理论分析的正确性和所得格式的精确性与有效性.考虑如下的二维变系数非齐次抛物型方程初边值问题其中:0<c1≤a(x,y)≤c2;Γ为Ω的边界.设τ=Δt=T/N为时间步长,h=Δx=Δy=1/M为空间步长,N、M均为正整数.表示在节点(jh, kh,nτ)处的网函数值,微分方程问题(1)—(3)的解函数为u(x,y,t),并记u(jh,kh,nτ)=u(j,k,n),由Taylor展开式在(j,k,n+1)处考虑方程(1),根据粘结系数法[7]有其中:分别为关于x、y的二阶中心差分算子;η1、η2、η3均为待定参数.适当选取这些参数使得式(9)—式(11)成立.由于即式(9)—式(11)的3个式子均成立,由式(8)、(14)、(16)和(17)得:将所得的等式左端加上右端加上,经过简单的计算可知,所加项的误差阶均不超过O(τ2+h4),则可得在(18)中略去小量项,并令,则可建立如下差分格式:易知差分格式(19)—(21)的截断误差为对格式(19)—(21)引进中间过渡量,可以得到问题(1)—(3)的ADI格式:当第n层上的值已知时,对于固定的k,(22)是关于的三对角线性方程组,可用追赶法求解.在解方程组(22)时,需要边值和(0≤k≤M-1),这可从方程(23)得到,在方程(23)中,分别令j=0和j=M,有当过渡层已求得时,对于固定的j,式(23)是关于的三对角线性方程组,也可用追赶法求解.解方程组(23)所需要的边界值和(0≤k≤M-1)可直接从已知的边界值条件得到,即:为对格式(22)—(23)进行稳定性分析,可将过渡层变量消去,即得格式(19),故两格式有相同的截断误差和稳定性质.因为差分格式关于初值稳定一定可以推出关于右端稳定,故只对齐次格式进行讨论,根据稳定性分析的Fourier方法[21],令将上式代入格式(19)式的齐次形式中,经计算整理,并利用关系式这里,可得格式(19)的传播因子为由s1的取值范围可知,所以|G(s1,s2)|根据Von neumann条件知格式(19)也即格式(22)—(23)恒稳定.综合上节与本节的论述,并根据Lax的稳定性与收敛性等价原理可得:定理本文的ADI格式恒稳定且以o(τ2+h4)的收敛阶收敛.在区域D:{0≤x,y≤1,t≥0}上对初值边问题用本文格式求数值解并与精确解u(x,y,t)=e-2tsin(x+y),相比较,取h=1/10,τ=rh2=r/100,r=1/2计算到n=200时的结果如表1所示.由表1结果看出,对所取的r,本文格式解与精确解均有很好的吻合,这与理论分析完全一致.本文构造了二维变系数抛物型方程的ADI格式,数值例子表明,它具有易于计算(能三对角追赶法求解)、精度高(较现有格式)、无条件稳定等特点.本文构造出的格式可以推广到任意维变系数抛物型方程,对当前关于变系数抛物型方程差分格式的构造与研究有重要的理论和实践意义.【相关文献】[1]DAI Weizhong,NASSAR Raja.A Second-order ADI scheme for three-dimensional parabolic differential equations[J].Numer Math,1964(6):428-453.[2]DAI Weizhong,NASSAR pact ADI method for solving parabolic differential equations[J].Numer Methods Partial Differential Eq,2002(18):129-142.[3] 孙志忠,李雪玲.反应扩散方程的紧交替方向差分格式[J].计算数学,2005,19(2):209-224.[4] 孙志忠,李雪玲.二维变系数反应扩散方程的紧交替方向差分格式[J].高等学校计算数学学报,2006,28(1):83-95.[5] 马菊意,杨辉.二维抛物型方程的一个高精度PC格式[J].工程数学学报,2008,25(2):373-376.[6]余德浩,汤华中.微分方程数值解法[M].北京:科学出版社,2003:102-174.[7] 戴嘉尊,邱建贤.微分方程数值解法[M].南京:东南大学出版社,2003:91-93[8]李立康,於崇华,朱政华.微分方程数值方法[M].上海:复旦大学出版社,1999:208-213.。
%%%%% 真解u=sin(pi*x)*sin(pi*y)*sin(t)
%%%%% 方程diff(u,t)-div(a(x,y)*grad(u))=f
%%%%% a(x,y)=x^2+y^2+1
%%%%%
f=sin(pi*x)*sin(pi*y)*cos(t)-2*pi^2*sin(pi*x)*sin(pi*y)*sin(t)*(x^2+y^2+1)+2*pi*sin(t)*(x*cos( pi*x)*sin(pi*y)+y*sin(pi*x)*cos(pi*y))
%clear all
% clc
%%%%finite element code for parabolic equation with constant coefficient
%%%mesh%%
node=[0,0;1,0;1,1;0,1];
elem=[2,3,1;4,1,3];
T=1;
bdEdge=setboundary(node,elem,’Dirichlet’);
n=input(‘Please input initial mesh:’);
M=input(‘M=’);
for i=1:n
[node,elem,bdEdge]=uniformrefine(node,elem,bdEdge);
end
N=size(node,1);
NT=size(elem,1);
S=1/NT;
r=1/M;
A=zeros(N,N);;
u=zeros(N,M+1);
u1=zeros(N,1);
f=inline(‘sin(pi*xx(1,1))*sin(pi*xx(1,2))*cos(t)-2*pi^2*sin(pi*xx(1,1))*sin(pi*xx(1,2))*sin(t)*(x x(1,1)^2+xx(1,2)^2+1)+2*pi*sin(t)*(xx(1,1)*cos(pi*xx(1,1))*sin(pi*xx(1,2))+y*sin(pi*xx(1,1))* cos(pi*xx(1,2));
a=inline(‘xx(1,1)^2+xx(1,2)^2+1’,’xx’);
[lambda,weight]=quadpts(5);
p=node’;
T=elem’;
totalEdge=[elem(:,[2,3]);elem(:,[3,1]);elem(:,[1,2])];
isBdEdge=reshape(bdEdge,3*NT,1);
Dirichlet=totalEdge(isBdEdge==1),:);
isBdNode=false(N,1);
isBdNode(Dirichlet)=true;
bdNode=find(isBdNode);
freeNode=find(~isBdNode);
for j=2:M+1
for i=1:NT
F=zeros(N,1);
F_ele=zeros(1,3);
T_ele=elem(i,:);
for m=1:7
xx(m,1)=(p(1,t(1,i))-p(1,t(3,i)))*lambda(m,1)+(p(1,t(2,i))-p(1,t(3,i))*lambda(m,2)+p(1,t(3,i));
xx(m,1)=(p(2,t(1,i))-p(2,t(3,i)))*lambda(m,1)+(p(2,t(2,i))-p(2,t(3,i))*lambda(m,2)+p(2,t(3,i));
end
for i=1:3
for k=1:7
F_ele(i)=F_ele(i)+S*weight(k)*lambda(k,i)*f(xx(k,:),(j-1)*r);
end
end
z=0;
for k=1:7
z=z+S*weight(k)*lambda(k)*a(xx(k,:));
end
x=node(T_ele,:);
[a,b]=Basis_coeff(x);
A_ele=[a,b]’*[a,b]/(4*s);
B_ele=zeros(3,3);
for i=1:3
for j=1:3
if i==j
B_ele(i,j)=1/12;
else
B_ele(i,j)=1/24;
end
end
end
A(T_ele,T_ele)=A(T_ele,T_ele)+S*B_ele+r*z*A_ele;
F(T_ele,1)=F(T_ele,1)+r*F_ele’+S*B_ele;
end
uj=zeros(N,1);
uj(freeNode)=A(freeNode,freeNode)\H(freeNode);
u(:,j)=uj;
end
showresult(node,elem,uj)
u_exact=zeros(N,M+1);
for j=1:M+1
u_exact(:,j)=inline(‘sin(pi*pxy(:,1)).*sin(pi*pxy(:,2)).*sin((j-1)*r)’,’pxy’);
end
L2_err=getL2error(node,elem,u_exact(:,j),u(:,j),5);
L’2_err=getL2error(node,elem,u_exact,u,5);
Du(:,j)=inline(‘[pi*cos(pi*pxy(:,1)).*sin(pi*pxy(:,2)).*sin((j-1)*r),pi*sin(pi*pxy(:,1)).*cos(pi*px y(:,2)).*sin((j-1)*r)]’,’pxy’);
H1_err=getH1error((node,elem,Dut(:,j),u(:,j),5);。