首页

偏微分方程数值解

偏微分方程数值解

课程设计

学院：专业：信息与计算科学学号：姓名：教师： xxxxxxxx 时间： 2011-12-3

一、问题提出

用Lax-Friedrich格式计算如下双曲型方程问题：

uu0,0x1,0t1.tx

u(x,0)cos(x),0x1.

u(0,t)u(1,t)cos(t),0t1.

其精确解为u(x,t)cos(xt)。二、网格剖分

为了用差分方法求解上述问题，将求解区域(x,t)|0x1,0t1作剖分。

]m等分，将时间[0,1]间作n等分，并记将空间区间[0,1作区

h1/m,

n1x/j,jh

。分别称h和为空间和时间步j,0mtkk,k,n0

长。用两簇平行直线xxj,0jm,ttk,0kn将分割成矩形网格。三、差分格式的建立 uu0………………………………（1） tx设G是x,t平面任一有界域，据Green公式：

uu

Gtx)dxdt(udtudx)

其中G。于是可将（1）式写成积分守恒形式：

(

(udtudx)0…………………………（2）



我们先从（2）式出发构造熟知的Lax格式设网格如下图所示

k+1

图一

取G为以A(j1,k)，B(j1,k1)，C(j1,k1)和D(j1,k)为顶点的矩形。

TABCDA为其边界，则

(udtudx)



DA

(u)dx(u)dx(u)dt(u)dt…………（3）

BC

AB

CD

右端第一个积分用梯形公式，第二个积分用中矩形公式，第三、四个积分用下矩形公式，则由（2）（3）式得到Lax-Friedrich格式

1k1

kkuk(uj1ukjj1)uuj1j10

2h

截断误差为

kk

Rkj(u)Lhuj[Lu]j

1k1k

kkkkuk(uujj1j1)uuuuj1j1jj

2htx



2k

uj

2t2

3k

h2ujO(h2),(0jm,0kn) 3

6x

所以Lax-Friedrich格式的截断误差的阶式O(h2) 令s/h：则可得差分格式为

sk1ksk1k

uku(uu)uj1,(0jm,0kn) jj1j1j1

222

u0jcos(xj)(0jm)

kk

u0cos(tk),umcos(tk),(0kn)

四、差分格式的求解

sk1ksk1ku(uu)uj1,(0jm,0kn) 差分格式ukjj1j1j1

222

写成如下矩阵形式：

1s220000

01s22000

1s220000



00

00

0

kk1

u0u1k10u1ku2



k1 k

uu

1sm2m1uk022m1

uk00m0

1s

22000



1s0220000

则需要通过对k时间层进行矩阵作用求出k+1时间层。

对上面的矩阵形式我通过C++编出如附录的程序求出数值解、真实解和误差。

对程序所得的结果进行分析如下表：

因为网格分的比较细，计算结果比较多，所以表1、表2、表3和表4分别给出不同网格得到的部分数值结果，数值解很好的逼近精确解。

表1

h1/90,1/100（s0.9)

表2

h1/900,1/1000(s0.9）

表3

h1/50,1/100（s0.5)

表4

h1/500,1/1000（s0.5)

当h(s1)时，数值解等于真实解，没有任何误差。当s取不同值时，误差曲线如图一、图二：

图一

图二

当s0.5时，网格疏密不同，误差曲线如下图

图三

当s0.9时，网格疏密不同，误差曲线如下图

图四

由图一和图二可以得出，当s越接近1（即和h接近或m和n）时， Lax-Friedrich格式的误差越小；由图三和图四可以得出，当s相同时，网格越密（即m和n越大），Lax-Friedrich格式的误差越小。

五、稳定性分析：

用Fourier方法对Lax-Friedrich格式进行稳定性分析：

k

令ukjve

ixj

代入上述的差分格式：

vk1e

six1ixsixix

vk[ej1(ej1ej1)ej1]

2221s

vk1vk[(eiheih)(eiheih)]vk(coshissinh)

22

ixj

G(s,h)coshissinh

G(s,h)cos2hs2sin2h1(s21)sin2h

2

s1，即



h

1时，此格式恒稳定。

附录

C++程序：

#include #include

#define p 50 //x方向的分段。

#define q 100 //y（时间）方向的分段。 #define pi 3.1415926

float A[p+1][q+1]; //定义二维数组A存放数值解。 float B[p+1][q+1]; //定义二维数组B存放真实解。

void main() {

float x[p+1],t[q+1],r,s,a,b; int i,j,k,n,m,w; r=1.0/p; s=1.0/q;

a=1.0/2-(1.0/2)*(s/r); //计算0.5-0.5*s。 b=1.0/2+(1.0/2)*(s/r); //计算0.5+0.5*s。

for(j=0;j

x[j]=j*r;

A[j][0]=B[j][0]=cos(pi*x[j]); }

for(k=1;k

t[k]=k*s;

A[0][k]=B[0][k]=cos(pi*t[k]); }

for(k=1;k

t[k]=k*s;

A[p][k]=B[p][k]=-cos(pi*t[k]); }

for(k=1;k

scanf("%d",&n); //输出第n时间层的数值解。 for(i=0;i

for(i=0;i

printf("%f ",B[i][n]); printf("\n");

for(i=0;i

printf("%f ",A[i][n]-B[i][n]);//输出数值解减去真实解即误差。 printf("\n");

scanf("%d,%d",&m,&w); //输出网格中任一点数值解、真实解和误差。 printf("%f %f %f",A[m][w],B[m][w],A[m][w]-B[m][w]); printf("\n"); }

Matlab程序：图一程序： t=0.1:0.1:1;

w1=[0.004251,0.016664,0.034592,0.052979,0.064844,0.067167,0.060746,0.047398,0.029052,0.007748];

w2=[0.000334,0.001320,0.002761,0.004362,0.005610,0.005533,0.004761,0.003521,0.001937,0.000163]; w3=[0,0,0,0,0,0,0,0,0,0];

plot(t,w1,'-bo',t,w2,'-k*',t,w3,'-rs')

图二程序：

t=0.1:0.1:1;

w1=[0.000454,0.001733,0.003580,0.005611,0.007148,0.007037,0.006015,0.004402,0.002359,0.000085];

w2=[0.000036,0.000138,0.000286,0.000448,0.000591,0.000560,0.000478,0.000348,0.000184,0.000002]; w3=[0,0,0,0,0,0,0,0,0,0];

plot(t,w1,'-bo',t,w2,'-k*',t,w3,'-rs'

图三程序：

t=0.1:0.1:1;

w1=[0.004251,0.016664,0.034592,0.052979,0.064844,0.067167,0.060746,0.047398,0.029052,0.007748];

w2=[0.000454,0.001733,0.003580,0.005611,0.007148,0.007037,0.006015,0.004402,0.002359,0.000085]; plot(t,w1,'-bd',t,w2,'-kp')

图四程序：

t=0.1:0.1:1;

w1=[0.000334,0.001320,0.002761,0.004362,0.005610,0.005533,0.004761,0.003521,0.001937,0.000163];

w2=[0.000036,0.000138,0.000286,0.000448,0.000591,0.000560,0.000478,0.000348,0.000184,0.000002]; plot(t,w1,'-bd',t,w2,'-kp')

偏微分方程数值解

课程设计

学院：专业：信息与计算科学学号：姓名：教师： xxxxxxxx 时间： 2011-12-3

一、问题提出

用Lax-Friedrich格式计算如下双曲型方程问题：

uu0,0x1,0t1.tx

u(x,0)cos(x),0x1.

u(0,t)u(1,t)cos(t),0t1.

其精确解为u(x,t)cos(xt)。二、网格剖分

为了用差分方法求解上述问题，将求解区域(x,t)|0x1,0t1作剖分。

]m等分，将时间[0,1]间作n等分，并记将空间区间[0,1作区

h1/m,

n1x/j,jh

。分别称h和为空间和时间步j,0mtkk,k,n0

长。用两簇平行直线xxj,0jm,ttk,0kn将分割成矩形网格。三、差分格式的建立 uu0………………………………（1） tx设G是x,t平面任一有界域，据Green公式：

uu

Gtx)dxdt(udtudx)

其中G。于是可将（1）式写成积分守恒形式：

(

(udtudx)0…………………………（2）



我们先从（2）式出发构造熟知的Lax格式设网格如下图所示

k+1

图一

取G为以A(j1,k)，B(j1,k1)，C(j1,k1)和D(j1,k)为顶点的矩形。

TABCDA为其边界，则

(udtudx)



DA

(u)dx(u)dx(u)dt(u)dt…………（3）

BC

AB

CD

右端第一个积分用梯形公式，第二个积分用中矩形公式，第三、四个积分用下矩形公式，则由（2）（3）式得到Lax-Friedrich格式

1k1

kkuk(uj1ukjj1)uuj1j10

2h

截断误差为

kk

Rkj(u)Lhuj[Lu]j

1k1k

kkkkuk(uujj1j1)uuuuj1j1jj

2htx



2k

uj

2t2

3k

h2ujO(h2),(0jm,0kn) 3

6x

所以Lax-Friedrich格式的截断误差的阶式O(h2) 令s/h：则可得差分格式为

sk1ksk1k

uku(uu)uj1,(0jm,0kn) jj1j1j1

222

u0jcos(xj)(0jm)

kk

u0cos(tk),umcos(tk),(0kn)

四、差分格式的求解

sk1ksk1ku(uu)uj1,(0jm,0kn) 差分格式ukjj1j1j1

222

写成如下矩阵形式：

1s220000

01s22000

1s220000



00

00

0

kk1

u0u1k10u1ku2



k1 k

uu

1sm2m1uk022m1

uk00m0

1s

22000



1s0220000

则需要通过对k时间层进行矩阵作用求出k+1时间层。

对上面的矩阵形式我通过C++编出如附录的程序求出数值解、真实解和误差。

对程序所得的结果进行分析如下表：

因为网格分的比较细，计算结果比较多，所以表1、表2、表3和表4分别给出不同网格得到的部分数值结果，数值解很好的逼近精确解。

表1

h1/90,1/100（s0.9)

表2

h1/900,1/1000(s0.9）

表3

h1/50,1/100（s0.5)

表4

h1/500,1/1000（s0.5)

当h(s1)时，数值解等于真实解，没有任何误差。当s取不同值时，误差曲线如图一、图二：

图一

图二

当s0.5时，网格疏密不同，误差曲线如下图

图三

当s0.9时，网格疏密不同，误差曲线如下图

图四

由图一和图二可以得出，当s越接近1（即和h接近或m和n）时， Lax-Friedrich格式的误差越小；由图三和图四可以得出，当s相同时，网格越密（即m和n越大），Lax-Friedrich格式的误差越小。

五、稳定性分析：

用Fourier方法对Lax-Friedrich格式进行稳定性分析：

k

令ukjve

ixj

代入上述的差分格式：

vk1e

six1ixsixix

vk[ej1(ej1ej1)ej1]

2221s

vk1vk[(eiheih)(eiheih)]vk(coshissinh)

22

ixj

G(s,h)coshissinh

G(s,h)cos2hs2sin2h1(s21)sin2h

2

s1，即



h

1时，此格式恒稳定。

附录

C++程序：

#include #include

#define p 50 //x方向的分段。

#define q 100 //y（时间）方向的分段。 #define pi 3.1415926

float A[p+1][q+1]; //定义二维数组A存放数值解。 float B[p+1][q+1]; //定义二维数组B存放真实解。

void main() {

float x[p+1],t[q+1],r,s,a,b; int i,j,k,n,m,w; r=1.0/p; s=1.0/q;

a=1.0/2-(1.0/2)*(s/r); //计算0.5-0.5*s。 b=1.0/2+(1.0/2)*(s/r); //计算0.5+0.5*s。

for(j=0;j

x[j]=j*r;

A[j][0]=B[j][0]=cos(pi*x[j]); }

for(k=1;k

t[k]=k*s;

A[0][k]=B[0][k]=cos(pi*t[k]); }

for(k=1;k

t[k]=k*s;

A[p][k]=B[p][k]=-cos(pi*t[k]); }

for(k=1;k

scanf("%d",&n); //输出第n时间层的数值解。 for(i=0;i

for(i=0;i

printf("%f ",B[i][n]); printf("\n");

for(i=0;i

printf("%f ",A[i][n]-B[i][n]);//输出数值解减去真实解即误差。 printf("\n");

scanf("%d,%d",&m,&w); //输出网格中任一点数值解、真实解和误差。 printf("%f %f %f",A[m][w],B[m][w],A[m][w]-B[m][w]); printf("\n"); }

Matlab程序：图一程序： t=0.1:0.1:1;

w1=[0.004251,0.016664,0.034592,0.052979,0.064844,0.067167,0.060746,0.047398,0.029052,0.007748];

w2=[0.000334,0.001320,0.002761,0.004362,0.005610,0.005533,0.004761,0.003521,0.001937,0.000163]; w3=[0,0,0,0,0,0,0,0,0,0];

plot(t,w1,'-bo',t,w2,'-k*',t,w3,'-rs')

图二程序：

t=0.1:0.1:1;

w1=[0.000454,0.001733,0.003580,0.005611,0.007148,0.007037,0.006015,0.004402,0.002359,0.000085];

w2=[0.000036,0.000138,0.000286,0.000448,0.000591,0.000560,0.000478,0.000348,0.000184,0.000002]; w3=[0,0,0,0,0,0,0,0,0,0];

plot(t,w1,'-bo',t,w2,'-k*',t,w3,'-rs'

图三程序：

t=0.1:0.1:1;

w1=[0.004251,0.016664,0.034592,0.052979,0.064844,0.067167,0.060746,0.047398,0.029052,0.007748];

w2=[0.000454,0.001733,0.003580,0.005611,0.007148,0.007037,0.006015,0.004402,0.002359,0.000085]; plot(t,w1,'-bd',t,w2,'-kp')

图四程序：

t=0.1:0.1:1;

w1=[0.000334,0.001320,0.002761,0.004362,0.005610,0.005533,0.004761,0.003521,0.001937,0.000163];

w2=[0.000036,0.000138,0.000286,0.000448,0.000591,0.000560,0.000478,0.000348,0.000184,0.000002]; plot(t,w1,'-bd',t,w2,'-kp')

相关文章

多种微分方程数值计算方法分析

2010年8月第4期城市勘测 Urban Geotechnical Investigation& Surveying Aug. 2010 No. 4 文章编号:1672-8262(2010)04-117-03中图分类号:O ...查看

对流扩散方程引言

对流扩散方程的定解问题是指物质输运与分子扩散的物理过程和黏性流体流动的数学模型,它可以用来描述河流污染.大气污染.核污染中污染物质的分布,流体的流动和流体中热传导等众多物理现象.关于对流扩散方程的求解很也备受关注,因此寻找一种稳定实用的数值 ...查看

常微分方程常用数值解法

第一章绪论 1.1 引言常微分方程是现代数学的一个重要分支,是人们解决各种实际问题的有效工具.微分方程的理论和方法从17世纪末开始发展起来,很快成了研究自然现象的强有力工具,在17到18世纪,在力学.天文.科学技术.物理中,就已借助微分 ...查看

如何看待数值模拟的作用

如何看待数值模拟的作用地下工程由于其特殊原因造就了其特殊性,与结构工程.机械工程等相比,有两大特点: 1. 地下工程所依赖的地质环境是地下工程的力学特征的主要影响因素: 2. 地质环境又是地下工程的承载单元. 因此,要对各种工程现象和力学 ...查看

数值分析学习心得体会

数值分析学习感想一个学期的数值分析,在老师的带领下,让我对这门课程有了深刻的理解和感悟.这门课程是一个十分重视算法和原理的学科,同时它能够将人的思维引入数学思考的模式,在处理问题的时候,可以合理适当的提出方案和假设.他的内容贴近实际, ...查看

基于三阶Adams格式的求解声波方程的多步算法

China University of Mining & Technology-Beijing 创新项目论文一种基于三阶Adams 格式的求解声波方程的多步算法摘要一个准确.高效.低数值频散的正演算法能够提高反演精度.加快反 ...查看

数学专业参考书推荐

数学专业参考书整理推荐从数学分析开始讲起: 数学分析是数学系最重要的一门课,经常一个点就会引申出今后的一门课,并且是今后数学系大部分课程的基础.也是初学时比较难的一门课,这里的难主要是对数学分析思想和方法的不适应,其实随着课程的深入会一点 ...查看

基于Matlab的传染病动力学模型仿真平台

基于Matlab的传染病动力学模型仿真平台 Simulation Platform of Epidemic Dynamics Model Based on MATLAB 摘要:开发了基于Matlab的传染病动力学模型仿真平台,通过对传染病动 ...查看

有限差分法.边界元法和离散元法

有限差分法已经发展的一些近似数值分析方法中,最初常用的是有限差分法,它可以处理一些相当困难的问题.但对于几何形状复杂的边界条件,其解的精度受到限制,甚至发生困难.作为60年代最重要的科技成就之一的有单元法.在理论和工程应用上都_得到迅速发 ...查看

算法大全第15章常微分方程的解法

第十五章常微分方程的解法建立微分方程只是解决问题的第一步,通常需要求出方程的解来说明实际现象,并加以检验.如果能得到解析形式的解固然是便于分析和应用的,但是我们知道,只有线性常系数微分方程,并且自由项是某些特殊类型的函数时,才可以肯定得 ...查看

热门内容