有限差分法求解偏微分方程MATLAB教学教材

有限差分法求解偏微分方程M A T L A B

南京理工大学

课程考核论文

课程名称：高等数值分析

论文题目：有限差分法求解偏微分方程姓名：罗晨

学号： 115104000545

成绩：

有限差分法求解偏微分方程

一、主要内容

1.有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程：

22(,)()u u

f x t t x

αα??-=??其中为常数

具体求解的偏微分方程如下：

22001

(,0)sin()(0,)(1,)00

u u x t x u x x u t u t t π???-=≤≤??????

=???

==≥???

2.推导五种差分格式、截断误差并分析其稳定性；

3.编写MATLAB 程序实现五种差分格式对偏微分方程的求解及误差分析；

4.结论及完成本次实验报告的感想。

二、推导几种差分格式的过程：

有限差分法（finite-difference methods ）是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。

推导差分方程的过程中需要用到的泰勒展开公式如下：

()2

100000000()()()()()()()......()(())

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 k t t x x h

τ

++-=??

-=? (2-2) 2.1 古典显格式

2.1.1 古典显格式的推导

由泰勒展开公式将(,)u x t 对时间展开得

2,(,)(,)(

)()(())i i k i k k k u

u 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 u

u x t u x t t t o t t t

u

u x t o t

ττ+++?=+-+-??=+?+? (2-4)

得到对时间的一阶偏导数

1,(,)(,)()=()i k i k i k u x t u x t u

o t ττ+-?+? (2-5) 由泰勒展开公式将(,)u x t 对位置展开得

223,,21(,)(,)()()()()(())2!k i k i k i i k i i u u

u 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,1122

231,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 i

u u

u 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,,22

231,,21(,)(,)()()()2!1(,)(,)()()()

2!i k i k i k i k i k i k i k i k

u u

u x t u x t h h o h x x u u u x t u x t h h o h x x +-???=+?+?+????????=-?+?+???? (2-8) 得到对位置的二阶偏导数

22

11,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)写成

1112

2k k

k k k k

i i i i i i u u u u u f h ατ

++-??--+-=????

(2-11) ()11

12

2k k k

k k k i i i i i i u u u

u 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)

2

r h

τ

=其中，古典显格式的差分格式的截断误差是2()o h τ+。

2.1.2 古典显格式稳定性分析

古典显格式(2-13)写成矩阵形式为

()112k k k h h h u ra I raC u f τ+=-++????

(2-14)

12212,(,,......,,)k k k k k

h N N r u u u u u h τ

--=

=其中。

(1)(1)

01

010*********N N C -?-??????

??=????????L L M

M L 上面的C 矩阵的特征值是：2cos()

1,2,......,1C j h j N λπ==-

()12H ra I raC =-+

()()()2

12=122cos()

121cos()14sin 1,2,......,1

2

H j C ra ra ra ra j h ra j h j h

ra j N λλπππ=-+-+=--=-=- (2-15)

使()1H ρ≤，即

2

114sin 12j h

ra π-≤-≤ 1

02ra ≤≤

结论：当1

02

ra ≤≤时，所以古典显格式是稳定的。

2.2 古典隐格式

2.2.1 古典隐格式的推导 将1k t t -=代入式 (2-3)得

21,11(,)(,)(

)()(())j k j k j k k k k k u

u x t u x t t t o t t t

---?=+-+-? (2-16) 21,(,)(,)()()j k j k j k u

u x t u x t o t

ττ-?=-?+? (2-17)

得到对时间的一阶偏导数

1,(,)(,)()=()j k j k j k u x t u x t u

o t ττ

--?+? (2-18) 将式(2-9)、(2-18)原方程得到

11122

(,)(,)

(,)2(,)(,)(,)()j k j k j k j k j k j k u x t u x t u x t u x t u x t f x t o h h αττ

-+---+??-=++????

(2-19)

为了方便把(2-19)写成

1

112

2k k k k k

j j

j j j k j u u u u u f h ατ

-+-??--+-=??????

(2-20) ()1

1

12

2k k k

k k k

j j

j j j j u u u

u u f h τα

τ-+----+= (2-21)

最后得到古典隐格式的差分格式为

()1

11(12)k k k k k j j j j

j ra u r u u u f ατ-+-+-+=+ (2-22) 2

r h

τ

=其中，古典隐格式的差分格式的截断误差是2()o h τ+。

2.2.2 古典隐格式稳定性分析

将古典隐格式(2-22)写成矩阵形式如下

()12

12()k k k

h h h

ra I raC u u f r h τ

τ++-=+=

???? (2-23)

误差传播方程

()112k k h

h ra I raC v v ++-=???? (2-24) ()12,A ra I raC B I

=+-=

所以误差方程的系数矩阵为

()1

112H A ra I raC --==+-????

()1

1,2,......,1122cos H j j N ra ra j h

λπ=

=-+-

使()1H ρ≤，显然

()2

1

122cos()1

12(1cos())

114sin 2

H j ra ra j h ra j h j h ra λπππ=

+-=+-=

+

1H j λ≤

恒成立。

结论：对于0r ?>，即任意网格比下，古典隐格式是绝对稳定的。

2.3 Richardson 格式

2.3.1 Richardson 格式的推导 将11k k t t t t +-==和，代入式(2-3)得

2

1,1121,11(,)(,)()()(())(,)(,)()()(())

i k i k i k k k k k i k i k i k k k k k

u u x t u x t t t o t t t u u x t u x t t t o t t t +++---??=+-+-????

??=+-+-???

(2-25) 即

2

1,21,(,)(,)()()(,)(,)()()

i k i k i k i k i k i k

u u x t u x t o t u u x t u x t o t ττττ+-??=+?+??????=-?+???

(2-26) 由此得到可得

211,(,)(,)()()2i k i k i k u x t u x t u

o t ττ++-?=+? (2-27)

将式(2-9) 、(2-27)代入原方程得到下式

22

11112(,)(,)(,)2(,)(,)(,)()2i 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-28) 为了方便可以把式(2-28)写成

11112

22k k k k k k

i i i i i i u u u u u f h ατ+-+-??--+-=????

(2-29) 即

()111

12

2k k k

k k k i i i i i i u u u

u u f h

τα

τ+-+---

-+= (2-30)

最后得到Richardson 显格式的差分格式为

()1111222k k k k k k i i i i i i u r u u u u f ατ+-+-=-+++ (2-31)

2

r h

τ

=其中，古典显格式的差分格式的截断误差是22()o h τ+。

2.3.2 Richardson 稳定性分析

将Richardson 显格式(2-31)写成如下矩阵形式

()11222k k k k h h h h u r C I u u f ατ+-=-++ (2-32)

误差传播方程矩阵形式

()1122k k k h h h

k

k

h h

v r C I v v v v α+-?=-+??=?? (2-33) 再将上面的方程组写成矩阵形式

112(2)0k k k h

k ra C I I v w

w I v ++-????== ????

??? (2-34) 系数矩阵的特征值是

4cos()4110j ra j h ra π-??∧=????

22

8sin 102

j h

ra πλλ+-= (2-35) 解得特征值为

1,2λ=

{

}2

12,=4sin 12j h Max ra πλλ> (恒成立) (2-37) 结论：上式对任意的网比都恒成立，即Richardson 格式是绝对不稳定的。

4. Crank-Nicholson 格式

3.4.1 Crank-Nicholson 格式的推导 将12

k t t +=代入式(2-9)得

111222111222

2,2+1,1+1+1(,)(,)()()(())(,)(,)()()(())

j j k j k k k k k k j j k j k k k k k k u u x t u x t t t o t t t u u x t u x t t t o t t t +++++++??

=+-+-????

??=+-+-???

(2-40) 即

12122

,2+1,1

(,)(,)()()2(,)(,)()()2j j k j k k j j k j k k u u x t u x t o t u u x t u x t o t ττττ+++??=+?+????

??=-?+???

(2-41)

得到如下方程

()

()

12122,1

2,12(,)(,)()=()

2(,)(,)()=()j j k k j k j j k k j k u x t u x t u o t u x t u x t u o t

ττττ++++???

-????+?????????

?

???-??

??-+?????????

(2-42) 所以12

k t t +=处的一阶偏导数可以用下式表示：

()(

)1122

12

,1,122(,)(,)2(,)(,)11()()()

22(,)(,)()

j j k j j k k k j k j k j k j k u x t u x t u x t u x t u u o t t u x t u x t o τττττ

+++++??--??????+=-+???

?

???????

?

-=+ (2-43)

将k t t =,11j j x x x x +-==和代入式(2-6)可以得到式(2-9)； 同理1k t t +=,11j j x x x x +-==和代入式(2-6)可以得到

2111112,122(,)2(,)(,)()()j k j k j k j k u x t u x t u x t u

o h x h

+++-++-+?=+? (2-44)

所以12

k t t +=，j x 处的二阶偏导数用式(2-6)、(2-44)表示：

22,,12211111112

22

1()()2(,)2(,)(,)(,)2(,)(,)1()2j k j k j k j k j k j k j k j k u u x x u x t u x t u x t u x t u x t u x t o h h h ++++-++-????+??????-+-+??=

++????

(2-45)

所以12

k t t +=，j x 处的函数值可用下式表示：

11

(,)(,)2j k j k f x t f x t +??+?? (2-46) 原方程变为：

()2212,1,,,12211()()()()()222k k j k j k j k j k j j u u u u

f f o h t t x x α+++????????+-+=++???????????? (2-47) 将差分格式代入上式得：

111111112

2221(,)(,)(,)2(,)(,)(,)2(,)(,)21(,)(,)()2

j k j k j k j k j k j k j k j k j k j k u x t u x t u x t u x t u x t u x t u x t u x t h h f x t f x t o h αττ++++-++-+--+-+??

-+????

??=+++?? (2-48) 为了方便写成：

()11111111112222k k k k k k k k k k j j j j j j j j j j r u u u u u u u u f f ατ++++++-+-????---++-+=+??????

(2-49)

最后得到Crank-Nicholson 格式的差分格式为

()()()()()1111111111=1+222k k k k k k k k j j j j j j j j r r r u u u r u u u f f αααατ+++++-+-??+-

+-+++????

(2-50) 2

r h τ

=其中，Crank-Nicholson 格式的差分格式的截断误差是22()o h τ+。

3.4.1 Crank-Nicholson 稳定性分析

Crank-Nicholson 格式写成矩阵形式如下：

()111(1)=(1)+222k k k k h h h h r r r I C u r I C u f f αααατ++??????+--++ ? ?????????

(2-51) 误差传播方程是：

1(1)=(1)22k k h h r r r I C v r I C v αααα+?

???+-

-+ ? ????

? (2-52) 可以得到：

1

(1)(1)22r r H r I C r I C αααα-?

???=+--+ ? ??

??? (2-53)

2

2

(1)cos()(1)cos()

12sin 212sin 2

H j r ra j h r ra j h j h ra j h

ra απλαπππ-+=

+--=

+ (2-54)

使()1H ρ≤ 即

2

2

12sin 21112sin 2

j h

ra j h

ra ππ--≤

≤+ (2-55) 222

12sin 12sin 12sin 222

j h j h j h

ra ra ra πππ--≤-≤+ (2-56) 2222

12sin 12sin 22

12sin 12sin 22j h j h ra ra j h j h ra ra ππππ?-≤+????--≤-?? (2-57)

22

222sin 2sin 2

2

sin 1sin 2

2

j h j h

ra ra j h j h ra ra ππππ?-≤???

?-≤-?? (2-58) 上式恒成立。

结论：Crank-Nicholson 格式对任意网格比也是绝对稳定的。

5. Du Fort Frankle 格式（Richardson 格式的改进）

将111

()2

k k k i i i u u u +-=+代入式(2-31)并化简得到Du Fort Frankle ：

()11111122k k k k k k k i i i i i i i u r u u u u u f ατ++--+-=--+++ (2-59) ()1111(12)2(12)2k k k k k i i i i i ra u r u u ra u f ατ+-+-+=++-+ (2-60) 可以证明Du Fort Frankle 格式是绝对稳定的。因为此格式是Richardson 格式的改进格式，因此截断误差还是22()o h τ+。

3.6 总结

(1) 古典显格式

古典显格式的差分格式为

()111(12)k k k k k i i i i i u ra u r u u f ατ++-=-+++

截断误差：2()o h τ+。

稳定性：当1

02

ra ≤≤时，古典显格式稳定。 (2) 古典隐格式

古典隐格式的差分格式为

()1

11(12)k k k k k j j j j

j ra u r u u u f ατ-+-+-+=+ 截断误差：2()o h τ+。

稳定性：任意网格比古典隐格式绝对稳定。 (3) Richardson 显格式

Richardson 显格式的差分格式为

()1111222k k k k k k i i i i i i u r u u u u f ατ+-+-=-+++

截断误差：22()o h τ+。

稳定性：任意网格比Richardson 格式绝对不稳定。

(4) Crank-Nicholson 格式

Crank-Nicholson 格式的差分格式为

()()()()()1111111111=1+222k k k k k k k k j j j j j j j j r r r u u u r u u u f f αααατ+++++-+-??+-

+-+++????

截断误差：22()o h τ+。

稳定性：Crank-Nicholson 格式对任意网格比绝对稳定。

(5) Du Fort Frankle 格式

()1111(12)2(12)2k k k k k i i i i i ra u r u u ra u f ατ+-+-+=++-+ 截断误差：22()o h τ+。

稳定性：Du Fort Frankle 格式对任意网格比绝对稳定。

三、MATLAB 实现五种差分格式对偏微分方程的求解及误差分析

3.1 精确数值解

22001

(,0)sin()(0,)(1,)00

u u x t x u x x u t u t t π???-=≤≤?????

=??==≥???

上述偏微分方程的精确解是

2

(,)sin()

(01,0)t u x t e x x t ππ-=≤≤≥

区域取值范围：01,00.2x t ≤≤≤≤。

用MATLAB 对精确解进行编程画出三维图像精确解程序如下：

close all clc

[x,t]=meshgrid(0:0.01:1,0:0.001:0.2)

u=exp(-pi*pi*t).*sin(pi*x) mesh(x,t,u); surf(x,t,u);

xlabel('x'),ylabel('t'),zlabel('u'); title('精确数值解'); shading interp

图3-1 精确数值解的三维图

(a) 精确数值解X-Y平面(b) 精确数值解X-Z平面

(c) 精确数值解Y-Z平面

图3-2 精确数值解的三个平面图

3.2 五种差分格式MATLAB程序

3.2.1 古典显格式

close all

clc

T=0.2

X=1.0

M=41

N=11

u=zeros(M,N); %构造一个M行N列的矩阵用于存放时间t和变量x

ra=(T/(M-1))/((X/(N-1))^2); %网格比

fprintf('稳定性系数 S=ra 为:\n');

disp(ra); % 显示网格比

for i=2:N-1

xx=(i-1)*(X/(N-1));

u(1,i)=sin(pi*xx);

end; %即t=0时刻赋值，边界条件

for k=1:M

u(k,1)=0;

u(k,N)=0;

end; % x=0,x=1处的边界条件

for k=1:M-1 %矩阵是从y轴表示行k，x轴表示列i。由行开始

for i=2:N-1

u(k+1,i)=((1-2*ra)*u(k,i)+ra*u(k,i+1)+ra*u(k,i-1)); %此处为古典显格式

end

end

disp(u); % 显示差分法求得的值

[x,t]=meshgrid(1:N,1:M); %将区域划分成网格对每个点赋值再画图

surf(x,t,u);

xlabel('x'),ylabel('t'),zlabel('u');

title(' 古典显格式'); %此程序得到的是图3-3

图3-3古典显格式程序结果图（r=0.5）

图3-4精确数值解、古典显格式程序结果的Y-Z平面图（r=0.5）

图3-5古典显格式在取不同网格比时的误差传播结果图

图3-6 不同时间取值时精确解、与古典显格式的值对比图（网格比r=0.5）红线

表示精确解、蓝色线表示差分格式的解

图3-1、图3-3对比可以看出，精确解和古典显格式（网格比r=0.5）的图形是一致的。图3-4精确数值解、古典显格式的Y-Z平面图结果可以看出古典显格式的边界值和精确解一样。图3-5是r分别取0.245、0.5、0.72、1.125时的误差传播图像，边界位置网格数为5处的误差为0.01得到的，可以看出r小于等于0.5是稳定的；但是r大于0.5出现不稳定现象。很好的验证了古典显格式稳定性。

3.2.2 古典隐格式

close all

clc

T=0.2

X=1.0

M=41

N=21

ra=(T/(M-1))/((X/(N-1))^2); %网格比

fprintf('稳定性系数 S=ra 为:\n');

disp(ra); %显示网格比

u=zeros(M,N); %构造一个M行N列的矩阵用于存放时间t和变量x

for i=2:N-1

xx=(i-1)*(X/(N-1));

u(1,i)=sin(pi*xx); %t=0时刻的赋值，边界条件

end;

for k=1:M

u(k,1)=0;

u(k,N)=0;

end; % x=0,x=1处的边界条件

A=zeros(N-1,N-1); %隐格式的矩阵形式中的A矩阵赋值

A(1,1)=1+2*ra;

A(N-1,N-1)=1+2*ra;

A(1,2)=-ra;

A(N-1,N-2)=-ra;

for m=2:N-2

A(m,m-1)=-ra;

A(m,m)=1+2*ra;

A(m,m+1)=-ra;

end;

%以下是——追赶法求u值

d=zeros(N-1,1); %隐格式右边初始矩阵

n=length(d);

U=zeros(n);

L=eye(n);

y=zeros(n,1);

x=zeros(n,1);

for i=1:N-1

d(i,1)=sin(pi*(i-1)*(1.0/(N-1)));

end %隐格式右边矩阵赋值

%以下循环对矩阵A进行LU分解

U(1,1)=A(1,1);

for i=2:n

L(i,i-1)=A(i,i-1)/U(i-1,i-1);

U(i-1,i)=A(i-1,i);%U的上次对角线即为A的上次对角线

U(i,i)=A(i,i)-L(i,i-1)*U(i-1,i);

end

for k=1:M-1 %外层大循环是求时间网格2,3……，M的求解u %以下是追赶法的求解过程

%---------追的过程--------即Ly=d的求解y

y(1,1)=d(1,1);

for i=2:n

y(i,1)=d(i,1)-L(i,i-1)*y(i-1,1);

end

%------赶的过程---------即Ux=y的求解x

x(n,1)=y(n,1)/U(n,n);

for i=n-1:-1:1

x(i,1)=(y(i,1)-U(i,i+1)*x(i+1,1))/U(i,i);

end %追赶法结束

for i=1:n

u(k+1,i)=x(i,1)

end

d=zeros(N-1,1); %更新右边矩阵

d=x %每次外循环更换右边矩阵

end

for k=1:M

u(k,1)=0;

end;

disp(u); % 显示差分法求得的值

[x,t]=meshgrid(1:N,1:M); %将区域划分成网格对每个点赋值再画图surf(x,t,u);

xlabel('x'),ylabel('t'),zlabel('u');

title('古典隐格式'); %此程序得到图是3-7

图3-7古典隐格式程序结果图（r=2）

图3-8精确数值解、古典隐格式程序结果的Y-Z平面图（r=2）

(完整版)偏微分方程的MATLAB解法

引言 偏微分方程定解问题有着广泛的应用背景。人们用偏微分方程来描述、解释或者预见各种自然现象，并用于科学和工程技术的各个领域fll。然而，对于广大应用工作者来说，从偏微分方程模型出发，使用有限元法或有限差分法求解都要耗费很大的工作量，才能得到数值解。现在，MATLAB PDEToolbox已实现对于空间二维问题高速、准确的求解过程。 偏微分方程 如果一个微分方程中出现的未知函数只含一个自变量，这个方程叫做常微分方程，也简称微分方程；如果一个微分方程中出现多元函数的偏导数，或者说如果未知函数和几个变量有关，而且方程中出现未知函数对几个变量的导数，那么这种微分方程就是偏微分方程。 常用的方法有变分法和有限差分法。变分法是把定解问题转化成变分问题，再求变分问题的近似解；有限差分法是把定解问题转化成代数方程，然后用计算机进行计算；还有一种更有意义的模拟法，它用另一个物理的问题实验研究来代替所研究某个物理问题的定解。虽然物理现象本质不同，但是抽象地表示在数学上是同一个定解问题，如研究某个不规则形状的物体里的稳定温度分布问题，由于求解比较困难，可作相应的静电场或稳恒电流场实验研究，测定场中各处的电势，从而也解决了所研究的稳定温度场中的温度分布问题。 随着物理科学所研究的现象在广度和深度两方面的扩展，偏微分方程的应用范围更广泛。从数学自身的角度看，偏微分方程的求解促使数学在函数论、变分法、级数展开、常微分方程、代数、微分几何等各方面进行发展。从这个角度说，偏微分方程变成了数学的中心。

一、MATLAB方法简介及应用 1.1 MATLAB简介 MATLAB是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。 1.2 Matlab主要功能 数值分析 数值和符号计算 工程与科学绘图 控制系统的设计与仿真 数字图像处理 数字信号处理 通讯系统设计与仿真 财务与金融工程 1.3 优势特点 1) 高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来； 2) 具有完备的图形处理功能，实现计算结果和编程的可视化； 3) 友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握； 4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ，

有限差分法求解偏微分方程MATLAB教学教材

有限差分法求解偏微分方程M A T L A B

南京理工大学 课程考核论文 课程名称：高等数值分析 论文题目：有限差分法求解偏微分方程姓名：罗晨 学号： 115104000545 成绩： 有限差分法求解偏微分方程

一、主要内容 1.有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程： 22(,)()u u f x t t x αα??-=??其中为常数 具体求解的偏微分方程如下： 22001 (,0)sin()(0,)(1,)00 u u x t x u x x u t u t t π???-=≤≤?????? =??? ==≥??? 2.推导五种差分格式、截断误差并分析其稳定性； 3.编写MATLAB 程序实现五种差分格式对偏微分方程的求解及误差分析； 4.结论及完成本次实验报告的感想。 二、推导几种差分格式的过程： 有限差分法（finite-difference methods ）是一种数值方法通过有限个微分方程近似求导从而寻求微分方程的近似解。有限差分法的基本思想是把连续的定解区域用有限个离散点构成的网格来代替；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。 推导差分方程的过程中需要用到的泰勒展开公式如下： ()2 100000000()()()()()()()......()(()) 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) 求解区域的网格划分步长参数如下：

Matlab求解微分方程(组)及偏微分方程(组)

第四讲 Matlab 求解微分方程（组） 理论介绍：Matlab 求解微分方程（组）命令 求解实例：Matlab 求解微分方程（组）实例 实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少.另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）.这就要求我们必须研究微分方程（组）的解法：解析解法和数值解法. 一．相关函数、命令及简介 1.在Matlab 中，用大写字母D 表示导数，Dy 表示y 关于自变量的一阶导数，D2y 表示y 关于自变量的二阶导数，依此类推.函数dsolve 用来解决常微分方程（组）的求解问题，调用格式为： X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解. 注意，系统缺省的自变量为t 2.函数dsolve 求解的是常微分方程的精确解法，也称为常微分方程的符号解.但是，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰富的函数，我们将其统称为solver ，其一般格式为： [T,Y]=solver(odefun,tspan,y0) 说明：(1)solver 为命令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一. (2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t

【毕业设计(论文)】二维热传导方程有限差分法的MATLAB实现

第1章前言 1.1问题背景 在史策教授的《一维热传导方程有限差分法的MATLAB实现》和曹刚教授的《一维偏微分方程的基本解》中，对偏微分方程的解得MATLAB实现问题进行过研究，但只停留在一维中，而实际中二维和三维的应用更加广泛。诸如粒子扩散或神经细胞的动作电位。也可以作为某些金融现象的模型，诸如布莱克-斯科尔斯模型与Ornstein-uhlenbeck过程。热方程及其非线性的推广形式也被应用与影响分析。 在科学和技术发展过程中，科学的理论和科学的实验一直是两种重要的科学方法和手段。虽然这两种科学方法都有十分重要的作用，但是一些研究对象往往由于他们的特性（例如太大或太小，太快或太慢）不能精确的用理论描述或用实验手段来实现。自从计算机出现和发展以来，模拟那些不容易观察到的现象，得到实际应用所需要的数值结果，解释各种现象的规律和基本性质。 科学计算在各门自然科学和技术科学与工程科学中其越来越大的作用，在很多重要领域中成为不可缺少的重要工具。而科学与工程计算中最重要的内容就是求解科学研究和工程技术中出现的各种各样的偏微分方程或方程组。 解偏微分方程已经成为科学与工程计算的核心内容，包括一些大型的计算和很多已经成为常规的计算。为什么它在当代能发挥这样大的作用呢？第一是计算机本身有了很大的发展；第二是数值求解方程的计算法有了很大的发展，这两者对人们计算能力的发展都是十分重要的。 1.2问题现状 近三十年来，解偏微分方程的理论和方法有了很大的发展，而且在各个学科技术的领域中应用也愈来愈广泛，在我国，偏微分方程数值解法作为一门课程，不但在计算数学专业，而且也在其他理工科专业的研究生的大学生中开设。同时，求解热传导方程的数值算法也取得巨大进展，特别是有限差分法方面，此算法的特点是在内边界处设计不同于整体的格式，将全局的隐式计算化为局部的分段隐式计算。而且精度上更好。 目前，在欧美各国MATLAB的使用十分普及。在大学的数学、工程和科学系科，MATLAB

有限差分法求解偏微分方程MATLAB

南京理工大学 课程考核论文 课程名称：高等数值分析 论文题目：有限差分法求解偏微分方程姓名：罗晨 学号： 成绩： 有限差分法求解偏微分方程

一、主要内容 1.有限差分法求解偏微分方程，偏微分方程如一般形式的一维抛物线型方程： 22(,)()u u f x t t x αα??-=??其中为常数 具体求解的偏微分方程如下： 22001 (,0)sin()(0,)(1,)00 u 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) 求解区域的网格划分步长参数如下：

Matlab求解微分方程(组)及偏微分方程(组)

第四讲Matlab求解微分方程(组) 理论介绍:Matlab求解微分方程(组)命令 求解实例:Matlab求解微分方程(组)实例 实际应用问题通过数学建模所归纳得到得方程,绝大多数都就是微分方程,真正能得到代数方程得机会很少、另一方面,能够求解得微分方程也就是十分有限得,特别就是高阶方程与偏微分方程(组)、这就要求我们必须研究微分方程(组)得解法:解析解法与数值解法、 一.相关函数、命令及简介 1、在Matlab中,用大写字母D表示导数,Dy表示y关于自变量得一阶导数,D2y 表示y关于自变量得二阶导数,依此类推、函数dsolve用来解决常微分方程(组)得求解问题,调用格式为: X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解、 注意,系统缺省得自变量为t 2、函数dsolve求解得就是常微分方程得精确解法,也称为常微分方程得符号解、但就是,有大量得常微分方程虽然从理论上讲,其解就是存在得,但我们却无法求出其解析解,此时,我们需要寻求方程得数值解,在求常微分方程数值解方 面,MATLAB具有丰富得函数,我们将其统称为solver,其一般格式为: [T,Y]=solver(odefun,tspan,y0) 说明:(1)solver为命令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb、ode15i之一、 (2)odefun就是显示微分方程在积分区间tspan上从到用初始条件求解、 (3)如果要获得微分方程问题在其她指定时间点上得解,则令tspan(要求就是单调得)、 (4)因为没有一种算法可以有效得解决所有得ODE问题,为此,Matlab提供了多种求解器solver,对于不同得ODE问题,采用不同得solver、 表1 Matlab中文本文件读写函数

MatlabPDE工具箱有限元法求解偏微分方程

在科学技术各领域中，有很多问题都可以归结为偏微分方程问题。在物理专业的力学、热学、电学、光学、近代物理课程中都可遇见偏微分方程。 偏微分方程，再加上边界条件、初始条件构成的数学模型，只有在很特殊情况下才可求得解析解。随着计算机技术的发展，采用数值计算方法，可以得到其数值解。 偏微分方程基本形式 而以上的偏微分方程都能利用PDE工具箱求解。 PDE工具箱 PDE工具箱的使用步骤体现了有限元法求解问题的基本思路，包括如下基本步骤： 1) 建立几何模型 2) 定义边界条件 3) 定义PDE类型和PDE系数 4) 三角形网格划分

5) 有限元求解 6) 解的图形表达 以上步骤充分体现在PDE工具箱的菜单栏和工具栏顺序上，如下 具体实现如下。 打开工具箱 输入pdetool可以打开偏微分方程求解工具箱，如下 首先需要选择应用模式，工具箱根据实际问题的不同提供了很多应用模式，用户可以基于适

当的模式进行建模和分析。 在Options菜单的Application菜单项下可以做选择，如下 或者直接在工具栏上选择，如下 列表框中各应用模式的意义为： ① Generic Scalar：一般标量模式（为默认选项）。 ② Generic System：一般系统模式。 ③ Structural Mech.，Plane Stress：结构力学平面应力。

④ Structural Mech.，Plane Strain：结构力学平面应变。 ⑤ Electrostatics：静电学。 ⑥ Magnetostatics：电磁学。 ⑦ Ac Power Electromagnetics：交流电电磁学。 ⑧ Conductive Media DC：直流导电介质。 ⑨ Heat Tranfer：热传导。 ⑩ Diffusion：扩散。 可以根据自己的具体问题做相应的选择，这里要求解偏微分方程，故使用默认值。此外，对于其他具体的工程应用模式，此工具箱已经发展到了Comsol Multiphysics软件，它提供了更强大的建模、求解功能。 另外，可以在菜单Options下做一些全局的设置，如下 l Grid：显示网格 l Grid Spacing…：控制网格的显示位置 l Snap：建模时捕捉网格节点，建模时可以打开 l Axes Limits…：设置坐标系范围 l Axes Equal：同Matlab的命令axes equal命令

一维导热方程 有限差分法 matlab实现

第五次作业（前三题写在作业纸上） 一、用有限差分方法求解一维非定常热传导方程，初始条件和边界条件见说明.pdf 文件，热扩散系数α=const ， 22T T t x α??=?? 1. 用Tylaor 展开法推导出FTCS 格式的差分方程 2. 讨论该方程的相容性和稳定性，并说明稳定性要求对求解差分方程的影响。 3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。 4. 编写M 文件求解上述方程，并用适当的文字对程序做出说明。（部分由网络搜索得到，添加，修改后得到。） function rechuandaopde %以下所用数据，除了t 的范围我根据题目要求取到了20000，其余均从pdf 中得来 a=0.00001;%a 的取值 xspan=[0 1];%x 的取值范围 tspan=[0 20000];%t 的取值范围 ngrid=[100 10];%分割的份数，前面的是t 轴的，后面的是x 轴的 f=@(x)0;%初值 g1=@(t)100;%边界条件一 g2=@(t)100;%边界条件二 [T,x,t]=pdesolution(a,f,g1,g2,xspan,tspan,ngrid);%计算所调用的函数 [x,t]=meshgrid(x,t); mesh(x,t,T);%画图，并且把坐标轴名称改为x ，t ，T xlabel('x') ylabel('t') zlabel('T') T%输出温度矩阵 dt=tspan(2)/ngrid(1);%t 步长 h3000=3000/dt;

h9000=9000/dt; h15000=15000/dt;%3000,9000,15000下，温度分别在T矩阵的哪些行T3000=T(h3000,:) T9000=T(h9000,:) T15000=T(h15000,:)%输出三个时间下的温度分布 %不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面 %稳定性讨论,傅里叶级数法 dx=xspan(2)/ngrid(2);%x步长 sta=4*a*dt/(dx^2)*(sin(pi/2))^2; if sta>0,sta<2 fprintf('\n%s\n','有稳定性') else fprintf('\n%s\n','没有稳定性') error end %真实值计算 [xe,te,Te]=truesolution(a,f,g1,g2,xspan,tspan,ngrid); [xe,te]=meshgrid(xe,te); mesh(xe,te,Te);%画图，并且把坐标轴名称改为xe，te，Te xlabel('xe') ylabel('te') zlabel('Te') Te%输出温度矩阵 %误差计算 jmax=1/dx+1;%网格点数 [rms]=wuchajisuan(T,Te,jmax) rms%输出误差

有限差分法的Matlab程序(椭圆型方程)

有限差分法的Matlab程序（椭圆型方程） function FD_PDE(fun,gun,a,b,c,d) % 用有限差分法求解矩形域上的Poisson方程 tol=10^(-6); % 误差界 N=1000; % 最大迭代次数 n=20; % x轴方向的网格数 m=20; % y轴方向的网格数 h=(b-a)/n; % x轴方向的步长 l=(d-c)/m; % y轴方向的步长 for i=1:n-1 x(i)=a+i*h; end % 定义网格点坐标 for j=1:m-1 y(j)=c+j*l; end % 定义网格点坐标 u=zeros(n-1,m-1); %对u赋初值 % 下面定义几个参数 r=h^2/l^2; s=2*(1+r); k=1; % 应用Gauss-Seidel法求解差分方程 while k<=N % 对靠近上边界的网格点进行处理 % 对左上角的网格点进行处理 z=(-h^2*fun(x(1),y(m-1))+gun(a,y(m-1))+r*gun(x(1),d)+r*u(1,m-2)+u(2,m-1))/s; norm=abs(z-u(1,m-1)); u(1,m-1)=z; % 对靠近上边界的除第一点和最后点外网格点进行处理 for i=2:n-2 z=(-h^2*fun(x(i),y(m-1))+r*gun(x(i),d)+r*u(i,m-2)+u(i+1,m-1)+u(i-1,m-1))/s; if abs(u(i,m-1)-z)>norm; norm=abs(u(i,m-1)-z); end u(i,m-1)=z; end % 对右上角的网格点进行处理 z=(-h^2*fun(x(n-1),y(m-1))+gun(b,y(m-1))+r*gun(x(n-1),d)+r*u(n-1,m-2)+u(n-2,m-1))/s; if abs(u(n-1,m-1)-z)>norm norm=abs(u(n-1,m-1)-z); end u(n-1,m-1)=z; % 对不靠近上下边界的网格点进行处理 for j=m-2:-1:2 % 对靠近左边界的网格点进行处理

电磁场实验一_有限差分法的matlab实现

电磁场与电磁波实验报告 实验项目：_______有限差分法__ ____ 班级：_____ __12电子2 ____ __ 实验日期：__2014年12月23日 姓名：___ _ __陈奋裕 __ __ 学号：___ ___1215106003 _____ 组员姓名：___ _ __ __ __ 组员学号：___ ___ _____ 指导教师：_ ____张海 ______

一、实验目的及要求 1、学习有限差分法的原理与计算步骤； 2、学习用有限差分法解静电场中简单的二维静电场边值问题； 3、学习用Matlab 语言描述电磁场与电磁波中内容，用matlab 求解问题并用图形表示出了，学习matlab 语言在电磁波与电磁场中的编程思路。 二、实验内容 理论学习：学习静电场中边值问题的数值法中的优先差分法的求解知识； 实践学习：学习用matlab 语言编写有限差分法计算二维静电场边值问题； 三、实验仪器或软件 电脑（WIN7）、Matlab7.11 四、实验原理 基本思想是把连续的定解区域用有限个离散点构成的网格来代替， 这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似， 积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组 ， 解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。 简单迭代法： 这一方法的求解过程是，先对场域内的节点赋予迭代初值(0),i j ?，这里上标(0)表示0次 （初始）近似值。然后按Laplace 方程 (k 1)(k)(k)(k)(k),1,,11,,11 []4 i j i j i j i j i j ?????+--++=+++(i,j=1,2,…) 进行反复迭代（k=0，1,2,…）。若当第N 次迭代以后，所有的内节点的相邻两次迭代值之间的最大误差不超过允许范围，即 (N)(N-1) ,,max|-|

Matlab求解微分方程(组)及偏微分方程(组)

第四讲 Matlab 求解微分方程（组） 理论介绍：Matlab 求解微分方程（组）命令 求解实例：Matlab 求解微分方程（组）实例 实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少.另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）.这就要求我们必须研究微分方程（组）的解法：解析解法和数值解法. 一．相关函数、命令及简介 1.在Matlab 中，用大写字母D 表示导数，Dy 表示y 关于自变量的一阶导数，D2y 表示y 关于自变量的二阶导数，依此类推.函数dsolve 用来解决常微分方程（组）的求解问题，调用格式为： X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解. 注意，系统缺省的自变量为t 2.函数dsolve 求解的是常微分方程的精确解法，也称为常微分方程的符号解.但是，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰富的函数，我们将其统称为solver ，其一般格式为： [T,Y]=solver(odefun,tspan,y0) 说明：(1)solver 为命令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一. (2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t 用初始条件0y 求解. (3)如果要获得微分方程问题在其他指定时间点012,,, ,f t t t t 上的解，则令 tspan 012[,,,]f t t t t =(要求是单调的). (4)因为没有一种算法可以有效的解决所有的ODE 问题，为此，Matlab 提供

《偏微分方程概述及运用matlab求解偏微分方程常见问题》要点

北京航空航天大学 偏微分方程概述及运用matlab求解微分方 程求解常见问题 姓名徐敏 学号57000211 班级380911班 2011年6月

偏微分方程概述及运用matlab求解偏微分 方程常见问题 徐敏 摘要偏微分方程简介，matlab偏微分方程工具箱应用简介，用这个工具箱解方程的过程是：确定待解的偏微分方程；确定边界条件；确定方程所在域的几何形状；划分有限元；解方程 关键词MATLAB 偏微分方程程序 如果一个微分方程中出现的未知函数只含有一个自变量，这个方程叫做常微分方程，也简称微分方程：如果一个微分方程中出现多元函数的偏导数，或者说如果未知函数和几个变量有关，而且方程中出现未知函数对几个变量的导数，那么这种微分方程就是偏微分方程。 一，偏微分方程概述 偏微分方程是反映有关的未知变量关于时间的导数和关于空间变量的导数之间制约关系的等式。许多领域中的数学模型都可以用偏微分方程来描述，很多重要的物理、力学等学科的基本方程本身就是偏微分方程。早在微积分理论刚形成后不久，人们就开始用偏微分方程来描述、解释或预见各种自然现象，并将所得到的研究方法和研究成果运用于各门科学和工程技术中，不断地取得了显著的成效，显示了偏微分方程对于人类认识自然界基本规律的重要性。逐渐地，以物

理、力学等各门科学中的实际问题为背景的偏微分方程的研究成为传统应用数学中的一个最主要的内容，它直接联系着众多自然现象和实际问题，不断地提出和产生出需要解决的新课题和新方法，不断地促进着许多相关数学分支(如泛函分析、微分几何、计算数学等)的发展，并从它们之中引进许多有力的解决问题的工具。偏微分方程已经成为当代数学中的一个重要的组成部分，是纯粹数学的许多分支和自然科学及工程技术等领域之间的一座重要的桥梁。 在国外，对偏微分方程的应用发展是相当重视的。很多大学和研究单位都有应用偏微分方程的研究集体，并得到国家工业、科学部门及军方、航空航天等方面的大力资助。比如在国际上有重大影响的美国的Courant研究所、法国的信息与自动化国立研究所等都集中了相当多的偏微分方程的研究人员，并把数学模型、数学方法、应用软件及实际应用融为一体，在解决实际课题、推动学科发展及加速培养人才等方面都起了很大的作用。 在我国，偏微分方程的研究起步较晚。但解放后，在党和国家的大力号召和积极支持下，我国偏微分方程的研究工作发展比较迅速，涌现出一批在这一领域中做出杰出工作的数学家，如谷超豪院士、李大潜院士等，并在一些研究方向上达到了国际先进水平。但总体来说，偏微分方程的研究队伍的组织和水平、研究工作的广度和深度与世界先进水平相比还有很大的差距。因此，我们必须继续努力，大力加强应用偏微分方程的研究，逐步缩小与世界先进水平的差距 二，偏微分方程的内容

Matlab解微分方程(ODE+PDE)

常微分方程： 1 ODE解算器简介(ode**) 2 微分方程转换 3 刚性/非刚性问题(Stiff/Nonstiff) 4 隐式微分方程(IDE) 5 微分代数方程(DAE) 6 延迟微分方程(DDE) 7 边值问题(BVP) 偏微分方程(PDEs)Matlab解法 偏微分方程： 1 一般偏微分方程组(PDEs)的命令行求解 2 特殊偏微分方程(PDEs)的PDEtool求解 3 陆君安《偏微分方程的MATLAB解法 先来认识下常微分方程(ODE)初值问题解算器(solver) [T,Y,TE,YE,IE] = odesolver(odefun,tspan,y0,options) sxint = deval(sol,xint) Matlab中提供了以下解算器： 输入参数： odefun：微分方程的Matlab语言描述函数，必须是函数句柄或者字符串，必须写成Matlab

规范格式(也就是一阶显示微分方程组)，这个具体在后面讲解 tspan=[t0 tf]或者[t0,t1,…tf]：微分变量的范围，两者都是根据t0和tf的值自动选择步长，只是前者返回所有计算点的微分值，而后者只返回指定的点的微分值，一定要注意对于后者tspan必须严格单调，还有就是两者数据存储时使用的内存不同(明显前者多)，其它没有任何本质的区别 y0=[y(0),y’(0),y’’(0)…]：微分方程初值，依次输入所有状态变量的初值，什么是状态变量在后面有介绍 options：微分优化参数，是一个结构体，使用odeset可以设置具体参数，详细内容查看帮助 输出参数： T：时间列向量，也就是ode**计算微分方程的值的点 Y：二维数组，第i列表示第i个状态变量的值，行数与T一致 在求解ODE时，我们还会用到deval()函数，deval的作用就是通过结构体solution计算t 对应x值，和polyval之类的很相似！ 参数格式如下： sol：就是上次调用ode**函数得道的结构体解 xint：需要计算的点，可以是标量或者向量，但是必须在tspan范围内 该函数的好处就是如果我想知道t=t0时的y值，不需要重新使用ode计算，而直接使用上次计算的得道solution就可以 [教程] 微分方程转换为一阶显示微分方程组方法 好，上面我们把Matlab中的常微分方程(ODE)的解算器讲解的差不多了，下面我们就具体开始介绍如何使用上面的知识吧！ 现实总是残酷的，要得到就必须先付出，不可能所有的ODE一拿来就可以直接使用，因此，在使用ODE解算器之前，我们需要做的第一步，也是最重要的一步，借助状态变量将微分

Matlab偏微分方程求解方法

Matlab 偏微分方程求解方法 目录： §1 Function Summary on page 10-87 §2 Initial Value Problems on page 10-88 §3 PDE Solver on page 10-89 §4 Integrator Options on page 10-92 §5 Examples” on page 10-93 §1 Function Summary 1.1 PDE Solver” on page 10-87 1,2 PDE Helper Functi on” on page 10-87 1.3 PDE Solver This is the MATLAB PDE solver. PDE Helper Function §2 Initial Value Problems pdepe solves systems of parabolic and elliptic PDEs in one spatial variable x and time t, of the form )x u ,u ,t ,x (s ))x u ,u ,t ,x (f x (x x t u )x u ,u ,t ,x (c m m ??+????=????- (10-2) The PDEs hold for b x a ,t t t f 0≤≤≤≤.The interval [a, b] must be finite. m

can be 0, 1, or 2, corresponding to slab, cylindrical, or spherical symmetry,respectively. If m > 0, thena ≥0 must also hold. In Equation 10-2,)x /u ,u ,t ,x (f ?? is a flux term and )x /u ,u ,t ,x (s ?? is a source term. The flux term must depend on x /u ??. The coupling of the partial derivatives with respect to time is restricted to multiplication by a diagonal matrix )x /u ,u ,t ,x (c ??. The diagonal elements of this matrix are either identically zero or positive. An element that is identically zero corresponds to an elliptic equation and otherwise to a parabolic equation. There must be at least one parabolic equation. An element of c that corresponds to a parabolic equation can vanish at isolated values of x if they are mesh points.Discontinuities in c and/or s due to material interfaces are permitted provided that a mesh point is placed at each interface. At the initial time t = t0, for all x the solution components satisfy initial conditions of the form )x (u )t ,x (u 00= (10-3) At the boundary x = a or x = b, for all t the solution components satisfy a boundary condition of the form 0)x u ,u ,t ,x (f )t ,x (q )u ,t ,x (p =??+ (10-4) q(x, t) is a diagonal matrix with elements that are either identically zero or never zero. Note that the boundary conditions are expressed in terms of the f rather than partial derivative of u with respect to x-x /u ??. Also, of

matlab实现有限差分法计算电场强度(最新)

实验一：有限差分法研究静电场边值问题 实验报告人：年级和班级：学号： 1. 实验用软件工具： Matlab 2. 实验原理：电磁场课本P36-38 1)差分方程 2)差分方程组的解 简单迭代法 高斯-赛德尔迭代法 逐次超松弛法 3. 实验步骤： 1）简单迭代法 程序： hx=41;hy=21; v1=zeros(hy,hx); v1(hy,:)=zeros(1,hx); v1(1,:)=ones(1,hx)*100; v1(:,1)=zeros(hy,1); v1(:,hx)=zeros(hy,1); v1 v2=v1;maxt=1;t=0; k=0; while(maxt>1e-5) k=k+1; maxt=0; for i=2:hy-1 for j=2:hx-1 v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v1(i-1,j)+v1(i,j-1))/4; t=abs(v2(i,j)-v1(i,j)); if(t>maxt) maxt=t;end end end v1=v2; end v2 k clf subplot(1,2,1),mesh(v2) axis([0,41,0,21,0,100]) subplot(1,2,2),contour(v2,15) hold on

axis([-1,42,-1,25]) plot([1,1,hx,hx,1],[1,hy+1,hy+1,1,1],'r') text(hx/2,0.3,'0V','fontsize',11); text(hx/2-0.5,hy+0.5,'100V','fontsize',11); text(-0.5,hy/2,'0V','fontsize',11); text(hx+0.3,hy/2,'0V','fontsize',11); hold off 当W=1e-5, 迭代次数：1401次 2）高斯-赛德尔迭代法 程序： hx=41;hy=21; v1=ones(hy,hx); v1(hy,:)=zeros(1,hx); v1(1,:)=ones(1,hx)*100; v1(:,1)=zeros(hy,1); v1(:,hx)=zeros(hy,1); v2=v1;maxt=1;t=0; k=0; while(maxt>1e-5) k=k+1; maxt=0; for i=2:hy-1 for j=2:hx-1 v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1))/4; t=abs(v2(i,j)-v1(i,j)); if(t>maxt) maxt=t;end end end v1=v2; end v2 k clf subplot(1,2,1),mesh(v2) axis([0,41,0,21,0,100]) subplot(1,2,2),contour(v2,15) hold on axis([-1,42,-1,25]) plot([1,1,hx,hx,1],[1,hy+1,hy+1,1,1],'r') text(hx/2,0.3,'0V','fontsize',11); text(hx/2-0.5,hy+0.5,'100V','fontsize',11); text(-0.5,hy/2,'0V','fontsize',11); text(hx+0.3,hy/2,'0V','fontsize',11); hold off

有限差分法

电磁场与电磁波项目训练报告 求解金属槽的电位分布 班级：通信13-4 姓名： 学号： 指导教师：徐维 成绩： 电子与信息工程学院 信息与通信工程系

求解金属槽的电位分布 1.实验原理 利用有限差分法和matlab软件解决电位在金属槽中的分布。 有限差分法基本思想是把连续的定解区域用有限个离散点构成的网格来代替,这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似,积分用积分和来近似,于是原微分方程和定解条件就近似地代之以代数方程组,即有限差分方程组,解此方程组就可以得到原问题在离散点上的近似解.然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解.在采用数值计算方法求解偏微分方程时,若将每一处导数由有限差分近似公式替代,从而把求解偏微分方程的问题转换成求解代数方程的问题。 2.有限差分法 方程的定解问题就是在满足某些定解条件下求微分方程的解。在空间区域的边界上要满足的定解条件称为边值条件。如果问题与时间有关，在初始时刻所要满足的定解条件，称为初值条件。不含时间而只带边值条件的定解问题，称为边值问题。与时间有关而只带初值条件的定解问题，称为初值问题。同时带有两种定解条件的问题，称为初值边值混合问题。 定解问题往往不具有解析解，或者其解析解不易计算。所以要采用可行的数值解法。有限差分方法就是一种数值解法，它的基本思想是先把问题的定义域进行网格剖分，然后在网格点上，按适当的数值微分公式把定解问题中的微商换成差商，从而把原问题离散化为差分格式，进而求出数值解。此外，还要研究差分格式的解的存在性和唯一性、解的求法、解法的数值稳定性、差分格式的解与原定解问题的真解的误差估计、差分格式的解当网格大小趋于零时是否趋于真解（即收敛性），等等。 有限差分方法具有简单、灵活以及通用性强等特点，容易在计算机上实现。 2.1有限差分法原理

一维扩散方程的有限差分法matlab

一维扩散方程的有限差分法 ——计算物理实验作业七 陈万 ● 题目： 编程求解一维扩散方程的解 取1.0,1.0,1.0,10,0.1,0,1,1,0,1,1max 0222111======-=====τh D t a c b a c b a 。输出t=1,2,...,10时刻的x 和u(x),并与解析解u=exp(x+0.1t)作比较。 ● 主程序： % 一维扩散方程的有限差分法 clear,clc; %定义初始常量 a1 = 1; b1 = 1; c1 = 0; a2 = 1;b2 = -1; c2 = 0; a0 = 1.0; t_max = 10; D = 0.1; h = 0.1; tao = 0.1; %调用扩散方程子函数求解 u = diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2); ● 子程序1： function output = diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2) % 一维扩散方程的有限差分法，采用隐式六点差分格式(Crank-Nicolson) % a0: x 的最大值 % t:_max: t 的最大值 % h: 空间步长 % tao: 时间步长

% D：扩散系数 % a1,b1,c1是（x=0）边界条件的系数；a2,b2,c2是（x=a0）边界条件的系数x = 0:h:a0; n = length(x); t = 0:tao:t_max; k = length(t); P = tao * D/h^2; P1 = 1/P + 1; P2 = 1/P - 1; u = zeros(k,n); %初始条件 u(1,:) = exp(x); %求A矩阵的对角元素d d = zeros(1,n); d(1,1) = b1*P1+h*a1; d(2:(n-1),1) = 2*P1; d(n,1) = b2*P1+h*a2; %求A矩阵的对角元素下面一行元素e e = -ones(1,n-1); e(1,n-1) = -b2; %求A矩阵的对角元素上面一行元素f f = -ones(1,n-1); f(1,1) = -b1; R = zeros(k,n);%求R %追赶法求解

五点差分法(matlab)解椭圆型偏微分方程

用差分法解椭圆型偏微分方程 U(0,y)=si n(pi*y),U(2,y)=eA2si n( pi*y); 先自己去看一下关于五点差分法的理论书籍 Matlab 程序： unction [p e u x y k]=wudianchafenfa(h,m,n,kmax,ep) % g-s 迭代法解五点差分法问题 %kmax 为最大迭代次数 %m,n 为x,y 方向的网格数，例如(2-0 ) /0.01=200; %e 为误差，p 为精确解 syms temp ; u=zeros(n+1,m+1); x=0+(0:m)*h; y=0+(0:n)*h; for (i=1:n+1) u(i,1)=sin(pi*y(i)); u(i,m+1)=exp(1)*exp(1)*sin(pi*y(i)); end for (i=1:n) for ( j=1:m) f(i,j)=(pi*pi-1)*exp(x(j))*sin(pi*y(i)); end -(Uxx+Uyy)=(pi*pi-1)eAxsin(pi*y) 0

end t=zeros(n-1,m-1); for (k=1:kmax) for (i=2:n) for ( j=2:m) temp=h*h*f(i,j)/4+(u(i,j+1)+u(i,j-1)+u(i+1,j)+u(i- 1,j))/4; t(i,j)=(temp-u(i,j))*(temp-u(i,j)); u(i,j)=temp; end end t(i,j)=sqrt(t(i,j)); if (k>kmax) break ; end if (max(max(t))