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

理工大学

课程考核论文

课程名称：高等数值分析

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

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

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

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

得到对位置的二阶偏导数

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

1112

2k k

k k k k

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

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

(2-11)

()11

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)

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

010100

101010N N C -?-??????

??=??

???

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

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

()12H ra I raC =-+

()()()2

12=122cos()

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

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

使()1H ρ≤，即

114sin 12j h

ra π-≤-≤ 1

02ra ≤≤

结论：当1

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

112

2k k k k k

j j

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

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

(2-20) ()11

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

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

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

21,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)代入原方程得到下式

2211112(,)(,)(,)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

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)

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)

误差传播方程矩阵形式

()11

22k k k h h h

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 π-??∧=????

8sin 102

j h

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

1,2λ=

(2-36)

{

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

4. Crank-Nicholson 格式

3.4.1 Crank-Nicholson 格式的推导将12

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

1112221112222

,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

,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

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 处的函数值可用下式表示：

(,)(,)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

222

1(,)(,)(,)2(,)(,)(,)2(,)(,)21(,)(,)()2j 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)22r r H r I C r I C αααα-?

???=+--+ ? ??

??? (2-53)

(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 ρ≤ 即

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 2212sin 12sin 22

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

(2-57)

222sin 2sin 2

sin 1sin 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

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 π???-=≤≤?????

=??==≥???

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

(,)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平面

图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）

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

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

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

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

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求解微分方程(组)及偏微分方程(组)

第四讲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

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

一、主要内容 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求解微分方程求解常见问题姓名徐敏学号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)解椭圆型偏微分方程

用差分法解椭圆型偏微分方程 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))

《matlab求解偏微分方程常见问题》

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

微分方程在MATLAB中的实现

微分方程在MATLAB中的实现作者：吴建宏时间：2011.09.20 ********************************************* 作者简介：吴建宏，男，毕业于哈尔滨工业大学（威海），现就读于同济大学，攻读汽车电子方向，有两年的MATLAB实践经验，个人喜欢建模，编程和电控方面的。真诚愿意和各位志同道合的朋友一起探讨交流，一起搭建广阔的知识平台。 *********************************************

PART ONE 微分方程简介一、微分方程基本概念微分方程：未知函数以及它们的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。其中，未知函数的最高阶数称为微分方程的阶。偏微分方程：如果未知函数是多元函数，那么这个微分方程就是偏微分方程。常微分方程：如果未知函数是一元函数，那么这个微分方程就是常微分方程。二、微分方程的解法高等数学里讲过，微分方程的解法有分离变量法，常数变易法，特征值法……而在解决工程实际问题时候，利用MATLAB 求解微分方程的方法主要有符号解法和数值解法。其中，符号解法主要可以求解可以用特征值法求解的常系数微分方程和少数特殊的微分方程，有很大的局限性，而且运行时间较长，优点就是可以得到精确的数学表达式。大部分实际的微分方程不得不通过数值解法来求解，优点就是运行时间快，可以解决复杂的线性或者非线性微分方程。缺点就是解是近似值。 PART TWO 微分方程的符号解法（dsolve ）微分方程的符号解法主要是函数dsolve ，这个函数用起来很简单，最重要的是要知道神马时候能够用dsolve ，神马时候不能用，当运行出错的时候，该怎么处理。接下来进入正文。一、语法 r =dsolve('eq1,eq2,...','cond1,cond2,...','v') r =dsolve('eq1','eq2',...,'cond1','cond2',...,'v') 二、语法详解 1.eq1,eq2……用来代表看上去能够用高等数学介绍过的手段求解出来的微分方程，如果实在不太熟悉或者忘了，没关系，后面会提供出现错误时候的处理措施；v 代表自变量，默认值为时间t 。边界条件/初始值用cond1,cond2,...来表示。 2.在eq1,eq2……中，用D 来代表变量的微分，例如?y 用Dy 表示，而? ?y 用y D 2来表示。 3.初始值的边界条件用b a y =)(或者b a Dy =)(来表示。如果初始值的个数小于变量的个数，输出的结果中就会出现待定系数。 4.Dsolve 能够接受的最大输入初始条件是12个 5.如果是一个等式一个输出，输出的结果是非线性的符号表达式；如果是多个等式和多个输出，将会按照左边输出参数[y1,y2……]中的顺序输出相应的符号表达式；如果是多等式，单输出，返回一个结构体数组。出现错误时候怎么处理？如果dsolve 不能找到有限解，它会试图去找方程的隐式解，当找到隐式解得时候，就会返回警告标志；如果既找不到有限解，也找不到隐式解，就会返回empty 的错误警告"warning Warning:explicit solution could not be found"。在这两种情况下都只能通过使用数值解法来解决。

MATLAB精通科学计算-偏微分方程求解

一、Maple V 系统 Maple V是由Waterloo大学开发的数学系统软件，它不但具有精确的数值处理功能，而且具有无以伦比的符号计算功能。Maple V的符号计算能力还是MathCAD和MATLAB等软件的符号处理的核心。Maple提供了2000余种数学函数，涉及范围包括：普通数学、高等数学、线性代数、数论、离散数学、图形学。它还提供了一套内置的编程语言，用户可以开发自己的应用程序，而且Maple自身的2000多种函数，基本上是用此语言开发的。 Maple采用字符行输入方式，输入时需要按照规定的格式输入，虽然与一般常见的数学格式不同，但灵活方便，也很容易理解。输出则可以选择字符方式和图形方式，产生的图形结果可以很方便地剪贴到Windows应用程序内。二、MATLAB 系统 MATLAB原是矩阵实验室（Matrix Laboratory）在70年代用来提供Linpack和Eispack 软件包的接口程序，采用C语言编写。从80年代出现3.0的DOS版本，逐渐成为科技计算、视图交互系统和程序语言。MATLAB可以运行在十几个操作平台上，比较常见的有基于Wi ndows 9X/NT、OS/2、Macintosh、Sun、Unix、Linux等平台的系统。 MATLAB程序主要由主程序和各种工具包组成，其中主程序包含数百个内部核心函数，工具包则包括复杂系统仿真、信号处理工具包、系统识别工具包、优化工具包、神经网络工具包、控制系统工具包、μ分析和综合工具包、样条工具包、符号数学工具包、图像处理工具包、统计工具包等。而且5.x版本还包含一套几十个的PDF文件，从MATLAB的使用入门到其他专题应用均有详细的介绍。 MATLAB是数值计算的先锋，它以矩阵作为基本数据单位，在应用线性代数、数理统计、自动控制、数字信号处理、动态系统仿真方面已经成为首选工具，同时也是科研工作人员和大学生、研究生进行科学研究的得力工具。MATLAB在输入方面也很方便，可以使用内部的Editor或者其他任何字符处理器，同时它还可以与Word6.0/7.0结合在一起，在Word的页面里直接调用MATLAB的大部分功能，使Word具有特殊的计算能力。三、MathCAD 系统 MathCAD是美国Mathsoft公司推出的一个交互式的数学系统软件。从早期的DOS下的1.0和Windows下的4.0版本，到今日的8.0版本，功能也从简单的数值计算，直至引用Maple强大的符号计算能力，使得它发生了一个质的飞跃。 MathCAD是集文本编辑、数学计算、程序编辑和仿真于一体的软件。MathCAD7.0 Pro fessional（专业版）运行在Win9X/NT下，它的主要特点是输入格式与人们习惯的数学书写格式很近似，采用WYSWYG（所见所得）界面，特别适合一般无须进行复杂编程或要求比较特殊的计算。MathCAD 7.0 Professional 还带有一个程序编辑器，对于一般比较短小，或者要求计算速度比较低时，采用它也是可以的。这个程序编辑器的优点是语法特别简单。 MathCAD可以看作是一个功能强大的计算器，没有很复杂的规则；同时它也可以和W ord、Lotus、WPS2000等字处理软件很好地配合使用，可以把它当作一个出色的全屏幕数学公式编辑器。四、Mathematica 系统 Mathematica是由美国物理学家Stephen Wolfram领导的Wolfram Research开发的数学系统软件。它拥有强大的数值计算和符号计算能力，在这一方面与Maple类似，但它的

(整理)Matlab解微分方程.

第十六章偏微分方程的数值解法科学研究和工程技术中的许多问题可建立偏微分方程的数学模型。包含多个自变量的微分方程称为偏微分方程(partial differential equation)，简称PDE 。偏微分方程问题，其求解是十分困难的。除少数特殊情况外，绝大多数情况均难以求出精确解。因此，近似解法就显得更为重要。本章仅介绍求解各类典型偏微分方程定解问题的差分方法。 16.1 几类偏微分方程的定解问题一个偏微分方程的表示通常如下： (,,,,)x x x y y y x y A B C f x y Φ+Φ+Φ=ΦΦΦ (16.1.1) 式中，,,A B C 是常数，称为拟线性(quasilinear)数。通常，存在3种拟线性方程：双曲型(hyperbolic)方程：240B AC ->；抛物线型(parabolic)方程：240B AC -=；椭圆型(ellliptic)方程：240B AC -<。 16.1.2 双曲型方程最简单形式为一阶双曲型方程： 0u u a t x ??+=?? (16.1.2) 物理中常见的一维振动与波动问题可用二阶波动方程： 22222u u a t x ??=?? (16.1.3) 描述，它是双曲型方程的典型形式。方程的初值问题为：

222220 0,(,0)()()t u u a t x t x u x x u x x t ?ψ=???=>-∞<<+∞ ????? =?? ??=-∞<<+∞ ??? (16.1.4) 边界条件一般有三类，最简单的初边值问题为： 22222120 00,0(,0)(0,)(),(,)()0()t u u a t T x l t x u x l u t g t u l t g t t T u x x t ?ψ=???==<<<?? (16.1.8) 方程可以有两种不同类型的定解问题： (1) 初值问题： 2200,(,0)()u u a t x t x u x x x ????-=>-∞<<+∞? ????=-∞<<+∞ ? (16.1.6) (2) 初边值问题： 2212 00,0(,0)()0(0,)(),(,)()0u u a t T x l t x u x x x l u t g t u l t g t t T ????-=<<<

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

用差分法解椭圆型偏微分方程 -（Uxx+Uyy)=(pi*pi-1)e^xsin(pi*y) 0kmax) break; end if(max(max(t))

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

[原创]偏微分方程数值解法的MATLAB源码【更新完毕】说明：由于偏微分的程序都比较长，比其他的算法稍复杂一些，所以另开一贴，专门上传偏微分的程序谢谢大家的支持！其他的数值算法见： ..//Announce/Announce.asp?BoardID=209&id=8245004 1、古典显式格式求解抛物型偏微分方程（一维热传导方程） function [U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C) %古典显式格式求解抛物型偏微分方程 %[U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C) % %方程：u_t=C*u_xx 0 <= x <= uX,0 <= t <= uT %初值条件：u(x,0)=phi(x) %边值条件：u(0,t)=psi1(t), u(uX,t)=psi2(t) % %输出参数：U -解矩阵，第一行表示初值，第一列和最后一列表示边值，第二行表示第2层…… % x -空间变量 % t -时间变量 %输入参数：uX -空间变量x的取值上限 % uT -时间变量t的取值上限 % phi -初值条件，定义为内联函数 % psi1 -边值条件，定义为内联函数 % psi2 -边值条件，定义为内联函数 % M -沿x轴的等分区间数 % N -沿t轴的等分区间数 % C -系数，默认情况下C=1 % %应用举例： %uX=1;uT=0.2;M=15;N=100;C=1; %phi=inline('sin(pi*x)');psi1=inline('0');psi2=inline('0'); %[U x t]=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C); %设置参数C的默认值 if nargin==7 C=1; end %计算步长 dx=uX/M;%x的步长 dt=uT/N;%t的步长

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 L 上的解,则令tspan 012[,,,]f t t t t =L (要求就是单调的)、 (4)因为没有一种算法可以有效的解决所有的ODE 问题,为此,Matlab 提供了多种求解器solver,对于不同的ODE 问题,采用不同的solver 、