当前位置:文档之家› 数理方程基于matlab的数值解法

数理方程基于matlab的数值解法

数理方程基于matlab的数值解法
数理方程基于matlab的数值解法

数理方程数值解法与其在matlab软件上的实现张体强1026222 廖荣发1026226

[摘要]

数学物理方程的数值解在实际生活中越来越使用,首先基于偏微分数值解的思想上,通过matlab软件的功能,研究其数学物理方程的数值解,并通过对精确解和数值解进行对比,追究其数值解的可行性,在此,给出相关例子和程序代码,利于以后的再次研究和直接使用。

[关键字]

偏微分方程数值解matlab 数学物理方程的可视化

一:研究意义

在我们解数学物理方程,理论上求数学物理方程的定解有着多种解法,但是有许多定解问题却不能严格求解,只能用数值方法求出满足实际需要的近似解。而且实际问题往往很复杂,这时即便要解出精确解就很困难,有时甚至不可能,另一方面,在建立数学模型时,我们已作了很多近似,所以求出的精确解也知识推导出的数学问题的精确解,并非真正实际问题的精确解。因此,我们有必要研究近似解法,只要使所求得的近似解与精确解之间的误差在规定的范围内,则仍能满足实际的需要,有限差分法和有限元法是两种最常用的

求解数学物理方程的数值解法,而MATLAB 在这一方面具有超强的数学功能,可以用来求其解。

二:数值解法思想和步骤

2.1:网格剖分

为了用差分方法求解上述问题,将求解区域

{}(,)|01,01x t x t Ω=≤≤≤≤作剖分。将空间区间[0,1]作m 等分,将时

间[0,1]区间作n 等分,并记

1/,1/,,0,,0j k h m n x jh j m t k k n ττ===≤≤=≤≤。分别称h 和τ

为空间和

时间步长。用两簇平行直线,0,,0j k x x j m t t k n =≤≤=≤≤将Ω分割成矩形网格。 2.2:差分格式的建立

0u u t x

??-=??………………………………(1) 设G 是,x t 平面任一有界域,据Green 公式(参考数学物理方程第五章):

(

)()G

u u

dxdt udt udx t x

Γ??-=--???

? 其中G Γ=?。于是可将(1)式写成积分守恒形式:

()0udt udx Γ

--=? (2)

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

图1

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

T ABCD = A 为其边界,则 ()()()()()DA

BC

AB

CD

udt udx u dx u dx u dt u dt Γ

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

??? (3)

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

1

11111()202k k k k k j j j j j u u u u u h

τ+-+-+-+-+=

截断误差为

()[]k k k

j h j j R u L u Lu =-

1

11111()22k k k k k k k j j j j j j j

u u u u u u u h t x

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

23222

3

(),(0,0)26k

k

j

j u u h O h j m k n t x

ττ??=

-=+≤≤≤≤?? 所以Lax-Friedrich 格式的截断误差的阶式2()O h τ+ 令/s h τ=:则可得差分格式为

1111

11(),(0,0)222

k k k k

k j j j j j s s u u u u u j m k n +--++=-+++≤≤≤≤ 0cos()(0)j j u x j m π=≤≤

k+1

0cos(),cos(),(0)k k k m k u t u t k n ππ==-≤≤

2.3差分格式的求解

差分格式111111(),(0,0)2

2

2

k k k k k j j j j j s s u u u u u j m k n +--++=-+++≤

≤≤≤

写成如下矩阵形式:

1

01112

12111100

0002222110000022221100000022220000000000

00

0k k k k k

k m m k m k m s s u u s s

u u u u

s s u u +++---??-+ ?

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

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

? ? ????? ? ??

?

则需要通过对k 时间层进行矩阵作用求出k+1时间层。 对上面的矩阵形式我通过C++(或matlab )编出如附录的程序求出数值解、真实解和误差。 例1:如下方程

0,01,0 1.(,0)cos(),0 1.(0,)(1,)cos(),0 1.

u u

x t t x

u x x x u t u t t t ππ??-=≤≤≤≤??=≤≤=-=≤≤ 利用 matlab 的数值解法结果并填入表中关系对比如下:

表1

1/900,1/1000(0.9h s τ===)

从上面可以看出,数值解精度相当高的

三:matlab的在数学物理方程上简单的应用

Matlab是一个强大的计算工具,超强的计算能力和绘图能力,下面几个例题来说明matlab数学物理上的应用

例1:将函数1/(1-a)2在z=0 处展成幂级数。

解:>>syms a;

>>Taylor (1/(1-a)^2,0)

ans =1+2*a+3*a^2+4*a^3+5*a^4+6*a^5

例2:写出函数f(x)=1/(x2+p2 )(a>0)的Fourier 变换式。

解:>>syms x w;

>>syms a potitive

>>f=1/(x^2+p^2);

>>F=Fourier (f,x,w)

F=pi*(signum(0,Re(p),0)*cosh(p*w)-2*Heaviside(w)*sinh(p

*w)+sinh(p*w))/p

例2:已知函数f(x)=x3 e-x,试求取该函数的Laplace 变换,并对结果进行Laplace 反变换。

解:>>syms x w;

>>f=x^3*exp(-x);

>>F=laplace(f,x,w)

F=6/(w+1)^4

对得出的结果进行Laplace 反变换,从而有

>>ilaplace(F)

ans=x^3*exp(-x)

利用手工方法对函数进行Fourier 变换和Laplace 变换,计算起来繁琐、复杂,且容易出错,利用MATLAB 快速、准确。

四:matlab解数学物理方程

4.1:数值解法与精确解的可视化对比分析

以下面的问题为例子(课本原题)

根据上面的可建立方程如下:

根据分离变量和差分其图形结果如下:

图2 分离变量热传导

图3 差分热传导

其代码如下:

图2

x=0:0.1*pi:pi;

y=0:0.4:10;

[x,t]=meshgrid(x,y);

u=0;

m=length(j);%matlab可计算的最大数,相当于无穷

for i=0:m

u=u+8*(-1)^i/(pi*(2*i+1)^2)*(sin((2*i+1)/2*x).*exp(-(2*i+1)^2/ 4*t));

end;

surf(x,t,u);

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

title(' 分离变量法(无穷)');

disp(u);

图3

u=zeros(20,100); %t=1 x=pi 20行100列横坐标为x 纵坐标为t s=(1/100)/(pi/20)^2;

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

disp(s);

for i=1:20

u(i,1)=i/20*pi;;

end;

for j=1:100

u(1,j)=0;

end

for j=1:99

for i=2:19

u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j);

end

end

for j=1:100

u(20,j)=u(19,j);

end;

disp(u);

[x,t]=meshgrid(1:100,1:20);

surf(x,t,u);

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

title(' 有限差分法解');

从上面可以看出数值解法精度很高,图形基本完全一样的

4.2 :matlab实现数值解法

以下面的方程为例

基本步骤:

(1)区域的离散或子区域的划分;

(2)插值函数的选择;

(3)方程组的建立;

(4)方程组的求解。

a=input(' 请输入系数a 的值:');

l=input(' 请输入长度l 的值:');

M=input(' 请输入将区间[0,l]等分的个数M:');

ot=input(' 请输入时间增量ot 的值:');

n=input(' 请输入运行次数n 的值:');

ox=l/M;x0=zeros(M+1,1);

for ii=1:M

x0(ii+1)=ii*ox;

end

u=sin(pi*x0/l);%t=0 时u(x,t)的值

r=a^2*ot/(ox)^2;

for ii=1:n

%数据的输入

B=zeros(M-1,1);%存放系数矩阵主对角线元素

A=zeros (M-2,1);%存放系数矩阵主对角线元素下方次对角线的元素

C=zeros (M-2,1);%存放系数矩阵主对角线元素上方次对角线的元素

S=zeros(M-1,1);%存放右端的常数项

for ii=1:M-2

B(ii)=1+2*r;A(ii)=-r;C(ii)=-r;

S(ii)=u(ii+1,1);

end

B(M-1)=1+2*r;S(M-1)=u(M,1);u(1,2)=0;u(M+1, 2)=0;

S(1,1)=S(1,1)+r*u(1,2);S(M-1,1)=S(M-1,1)+r*u (M+1,2);

%追赶法

S(1)=S(1)/B(1);T=B(1);k=2;

while k~=M

B(k-1)=C(k-1)/T;

T=B(k)-A(k-1)*B(k-1);

S(k)=(S(k)-A(k-1)*S(k-1))/T;

k=k+1;

end

k=1;

while k~=M-1

S(M-1-k)=S(M-1-k)-B(M-1-k)*S(M-k);

k=k+1;

end

u(2:M,2)=S; %把结果放入矩阵u 中

u(:,1)=u(:,2);% 过河拆桥

end

%计算精确值,存放在u 的第二列

for x=0:M

u(x+1,2)=exp(-(pi*a/l)^2*n*ot)*sin(pi*x*ox/l); end

%计算最大相对误差

ez=zeros(M-1,1);

for ii=2:M

ez(ii-1)=abs(u(ii,1)-u(ii,2))/u(ii,2);

end

E=max(ez);

fprintf (' 最后时刻数值解与精确解分别为:\n');disp (u);

fprintf (' 差分法得到的结果与正确结果的最大相对误差为:');

disp([num2str(E*100) '%']);

%画二维图比较

plot(x0,u(:,1),'r:',x0,u(:,2),'b-');

legend(' 数值解',' 精确解')

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

title(' 最后时刻热传导问题数值解与精确解比较')。运行如下:

请输入系数a 的值:1

请输入长度l 的值:1

请输入将区间[0,l]等分的个数M:10 请输入时间增量ot 的值:0.001

请输入运行次数n 的值:100

最后时刻数值解与精确解分别为:

0 0

0.1167 0.1152

0.2219 0.2191

0.3054 0.3015

0.3591 0.3545

0.3775 0.3727

0.3591 0.3545

0.3054 0.3015

0.2219 0.2191

0.1167 0.1152

0 0.0000

差分法得到的结果与正确结果的最大相对误差为:1.2934%

图4 数值解与精确解对比

五: matlab偏微分工具箱介绍

上面我们用的是编程技术,和c++没有多大区别,但是在matlab中还提供了偏微分工具箱,可以模拟数理方程动态图,

这里简单介绍一些功能

5.1:偏微分方程工具箱的功能

偏微分方程工具箱(PDE Toolbox)提供了研究和求解空间二维偏微分方程问题的一个强大而又灵活实用的环境。PDE Toolbox的功能包括:

(1) 设置PDE (偏微分方程)定解问题,即设置二维定解区域、边界条件以及方程的形式和系数;

(2) 用有限元法(FEM) 求解PDE数值解;

(3) 解的可视化。

下面我简单介绍一下操作

启动matlab,在工作区输入pdetool,会出现一个界面,如下

图5 PDE Toolbox界面

系统立即产生偏微分方程工具箱(PDE Toolbox)的图形用户界面(Graphical User Interface,简记为GUI),即PDE解的图形环

境,这时就可以在它上面画出定解区域、设置方程和边界条件、作网格剖分、求解、作图等工作。

5.2:偏微分工具箱的求解范围

5.2.1:方程类型

PDE Toolbox求解的基本方程有椭圆型方程、抛物型方程、双曲型方程、特征值方程、椭圆型方程组以及非线性椭圆型方程。(参考这学期数学物理方程几大类型)另外,对于Poission方程还有一个矩形网格的快速求解器

5.2.2 如何使用FDE Toolbox

1:定解问题的设置

员简单的办法是在PDE Tool上直接使用图形用户界面(GUl)。设置定解问题包括三个步骤:

(1)Draw模式:使用CSG(几何结构实体模型)对话框画几何区域,包括矩形、圆、椭圆和多边形,也可以将它们组合使用。

(2)Boundary模式:在各个边界段上给出边界条件,

(3)PDE模式:确定方程的类型、系数c,a,f和d c。也能够在不同子区域上设置不同的系数(反映材料的性质)。

5.2.3 解PDE问题

用GUI解PDE问题主要经过下面两个过程(模式)

(1)Mesh模式;生成网格.自动控制网格参数。

(2)Solve模式:对于椭圆型方程还能求非线性和自适应解。对于抛物型和双曲型力程.设置初始边值条件后能求出给定t时刻的解。

对于特征值问题,能求出给定区间内的特征值;求解后可以加密网格再求解。

5.2.3 使用Toolbox求解非标准的问题

对于非标准的问题。可以用PDE Too1box的函数。或者用FEM(有限元法)求解更为复杂的问题。

5.2.4 计算结果的可视化

从GUI能够使用Plot模式实现可视化。可以使用Color, Height 和Vector等作图。对于抛物型和双曲型方程,还可以生成解的动画。这些操作通过命令行都很容易实现。

5.2.5 应用领域

在应用界面提供了丁如下应用领域

.结构力学——平面应力问题

.结构力学——平面应变问题

.静电场问题

.静磁场问题

.交流电磁场问题

.直流导体介质问题

.热传导问题

下面是一个我们模拟一个简单的二维波动方程,其中程序部分如下[p,e,t]=initmesh(‘squareg’);

x=p(1,:)’;y=p(2,:)’;

u0=atan(cos(pi/2*x));

ut0=3*sin(pi*x).*exp(sin(pi/2*y));

n=3l;

tlist=linspace(0,5,n);

uu=hyperbolic(u0,ut0,tlist,‘squareb3’,p,e,t,1,0,0,1);

delta=-1:0.1:1;

[uxy,tn,a2,a3]=tri2grid(p,t,uu(:,1),delta,delta);

gp=[tn;a2;a3];

umax=max(max(uu));

umin=min(min(uu));

newplot

M=moviein(n);

For I=1:n,

pdeplot(p,e,t,‘xydata’,uu(:,I),‘zdata’,uu(:,I), ‘mesh’,‘off’,‘xygrid’, ‘on’, ‘gridparam’,gp,‘colorbar’,‘off’,‘zstyle’, ‘continuous’);

axis([-1 1 –1 1 umin umax]);caxis([umin umax]);

M(:,I)=getframe;

End

Movie(M,10);

图6 matlab模拟二维波动方程

六:总结

在以后我们现实生活中很少出现非数值解,数值解法的思想很重要,在基于思想上我们可以通过计算机实现,matlab 在这一方面非常方便,可以通过纯的编程实现其解的值,也可以借助其工具箱解方程,可以模拟二维图,三维图和动态图,学习不但只是学,还要懂得和生活练习,去实践,而本文章通过课本上的例题进而分析数值解法与精确解的区别,从而认为数值解不但可行,而且很使用,在计算机实现数学物理方程的数值解时,只要对以上的程序稍作修改,即可得数值解,还可以通过工具箱模拟出来,为科学研究提供一个

更好的平台。

参考文献

[1]曹钢,王桂珍,任晓荣.一维热传导方程的基本解[J].山东轻

工业学院学报, 2005,19(4):76-80.

[2]万正苏,方春华,张再云.关于热传导方程有限差分区域分解并行算法精度的注记[J].湖南理工学院学报(自然科学版),2007 ,20(3):12-14.

[3]StephenJ.Chapman.MATLAB 编程[M].邢树军,郑碧波,译. 北京:科学出版社,2008.

[4]田兵.用MATLAB 解偏微分方程[J].阴山学刊,2006,20(4): 12-13.

[5]王飞,裴永祥.有限差分方法的MATLAB 编程[J].新疆师

范大学学报(自然科学版),2003,22(4):21-27.

[6]王宝红.热传导方程的可视化探讨[J].忻州师范学院学报, 2008,24(2):31-36.

[7]李先枝.热传导方程差分解法的最佳网格[J].河南大学学报(自然科学版),2004,34(3):16-18.

[8]赵德奎,刘勇.MATLAB 在有限差分数值计算中的应用[J]. 四川理工学院学报,2005,18(4):61-64.

[9]谢焕田,吴艳.拉普拉斯有限差分法的MATLAB 实现[J].

四川理工学院学报,2008,21(3):1-2.

[10] 赵刚,严尚安等主编,《数学物理方程》,中国科技出版社,2004年版,第8页

[11] 复旦大学数学系主编,《数学物理方程》,上海科学技术出版社,1960年版,第10页

[12] 闻新等编著,《MTLAB科学图形构建基础与应用》,科学出版社,2002年版,第107页

[13] Zakarauskas Rierre, Dosso Stanley E, Fawcett John A. Matched-field inversion for source location and optimal iquivalent bathymetry [J]. J

[14] 梁昆淼编,《数学物理方法》,人民教育出版社1978年

第二版,第256—257页,第390—392页

[15] 彭芳麟著,《数学物理方程的MATLAB解法与可视化》,

清华大学出版社,2004年版,第140页,第142页,第150页

[16] 李好,杨天春,王齐仁. 基于Matlab7. 0 PDE工具箱求解数学物理方程[ J ]. 电脑开发与应用, 2009, 22 (1) : 26 - 27, 48.

2012年5月

线性代数方程组数值解法及MATLAB实现综述

线性代数方程组数值解法及MATLAB实现综述廖淑芳20122090 数计学院12计算机科学与技术1班(职教本科)一、分析课题 随着科学技术的发展,提出了大量复杂的数值计算问题,在建立电子计算机成为数值计算的主要工具以后,它以数字计算机求解数学问题的理论和方法为研究对象。其数值计算中线性代数方程的求解问题就广泛应用于各种工程技术方面。因此在各种数据处理中,线性代数方程组的求解是最常见的问题之一。关于线性代数方程组的数值解法一般分为两大类:直接法和迭代法。 直接法就是经过有限步算术运算,可求的线性方程组精确解的方法(若计算过程没有舍入误差),但实际犹如舍入误差的存在和影响,这种方法也只能求得近似解,这类方法是解低阶稠密矩阵方程组级某些大型稀疏矩阵方程组的有效方法。直接法包括高斯消元法,矩阵三角分解法、追赶法、平方根法。 迭代法就是利用某种极限过程去逐步逼近线性方程组精确解的方法。迭代法具有需要计算机的存储单元少,程序设计简单,原始系数矩阵在计算过程始终不变等优点,但存在收敛性级收敛速度问题。迭代法是解大型稀疏矩阵方程组(尤其是微分方程离散后得到的大型方程组)的重要方法。迭代法包括Jacobi法SOR法、SSOR法等多种方法。 二、研究课题-线性代数方程组数值解法 一、直接法 1、Gauss消元法 通过一系列的加减消元运算,也就是代数中的加减消去法,以使A对角线以下的元素化为零,将方程组化为上三角矩阵;然后,再逐一回代求解出x向量。

1.1消元过程 1. 高斯消元法(加减消元):首先将A 化为上三角阵,再回代求解。 11121121222212n n n n nn n a a a b a a a b a a a b ?? ? ? ? ??? (1)(1)(1)(1)(1)11 121311(2)(2)(2)(2)222322 (3)(3)(3)3333()()000 00 n n n n n nn n a a a a b a a a b a a b a b ?? ? ? ? ? ? ?? ? 步骤如下: 第一步:1 11 1,2,,i a i i n a -? +=第行第行 11121121222212 n n n n nn n a a a b a a a b a a a b ?? ? ? ? ??? 111211(2)(2)(2)2222 (2)(2)(2)2 00n n n nn n a a a b a a b a a b ?? ? ? ? ??? 第二步:(2)2 (2)222,3, ,i a i i n a -?+=第行第行 111211(2)(2)(2)2222 (2)(2)(2)2 00n n n nn n a a a b a a b a a b ?? ? ? ? ??? 11121311(2)(2)(2)(2)222322 (3)(3)(3)33 33(3)(3)(3)3 00000 n n n n nn n a a a a b a a a b a a b a a b ?? ? ? ? ? ? ?? ? 类似的做下去,我们有: 第k 步:() ()k ,1, ,k ik k kk a i i k n a -?+=+第行第行。 n -1步以后,我们可以得到变换后的矩阵为: 11121311(2)(2)(2)(2)222322 (3)(3)(3)3333()()00000 n n n n n nn n a a a a b a a a b a a b a b ?? ? ? ? ? ? ?? ?

有限差分法求解偏微分方程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 软件求解微分方程 1.解析解 (1)一阶微分方程 求21y dx dy +=的通解:dsolve('Dy=1+y^2','x') 求y x dx dy -+=21的通解:dsolve('Dy=1+x^2-y','x') 求?????=+=1 )0(12y y dx dy 的特解:dsolve('Dy=1+y^2',’y(0)=1’,'x') (2)高阶微分方程 求解???-='==-+'+''. 2)2(,2)2(,0)(222πππy y y n x y x y x 其中,21=n ,命令为: dsolve('x^2*D2y+x*Dy+(x^2-0.5^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x') 求042=-+'-'''x y y y 的通解,命令为: dsolve('D3y-2*Dy+y-4*x=0','x') 输出为: ans=8+4*x+C1*exp(x)+C2*exp(-1/2*(5^(1/2)+1)*x)+C3*exp(1/2*(5^(1/2)-1)*x) (3)一阶微分方程组 求???+-='+='). (3)(4)(),(4)(3)(x g x f x g x g x f x f 的通解:[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','x') 输出为: f =exp(3*x)*(cos(4*x)*C1+sin(4*x)*C2) g =-exp(3*x)*(sin(4*x)*C1-cos(4*x)*C2) 若再加上初始条件1)0(,0)0(==g f ,则求特解: [f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0,g(0)=1','x') 输出为: f =exp(3*x)*sin(4*x) g =exp(3*x)*cos(4*x) 2.数值解 (1)一阶微分方程

(完整版)偏微分方程的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就是显示微分方程在积分区间tspan上从到用初始条件求解、 (3)如果要获得微分方程问题在其她指定时间点上得解,则令tspan(要求就是单调得)、 (4)因为没有一种算法可以有效得解决所有得ODE问题,为此,Matlab提供了多种求解器solver,对于不同得ODE问题,采用不同得solver、 表1 Matlab中文本文件读写函数

matlab数值计算(命令与示例)

MATLAB数值计算 MATLAB数值计算 (1) 1创建矩阵 (3) 1.1直接输入 (3) 1.2向量 (3) 1.2.1linspace:线性分布 (3) 1.2.2冒号法 (3) 1.3函数创建 (4) 1.3.1eye:单位矩阵 (4) 1.3.2rand:随机矩阵 (4)

1.3.3zeros:全0矩阵 (4) 1.3.4ones:全1矩阵 (5) 2矩阵运算 (5) 2.1加减 (5) 2.1.1[M×N]±[M×N] (5) 2.2乘 (6) 2.2.1[M×N]*a (6) 2.2.2[M×N]*[N×M] (6) 2.3乘方 (7) 2.3.1[M×M]^a (7) 2.3.2a^[M×M] (7) 2.4特殊运算 (8) 2.4.1求逆inv (8) 2.4.2行列式det (8) 2.4.3特征值eig (8) 2.4.4转置'和.' (9) 2.4.5变形reshape (10) 2.4.6翻转rot90,fliplr,flipud (11) 2.4.7抽取diag,tril,triu (12) 2.5数组运算 (12) 2.5.1乘 (12) [M×N].*[M×N] (12) 2.5.2除 (13) [M×N]./[M×N] (14) [M×N].\[M×N] (14) 2.5.3乘方 (14) [M×N].^[M×N] (15) a.^[M×N] (15) 2.6除法 (15) 2.6.1求解线性方程组 (15) 3多项式 (16) 3.1系数表示法poly (16) 3.2求根roots (16) 3.3乘法conv (16) 3.4除法deconv (17) 3.5求值polyval (17) 3.6微分polyder (18)

一维抛物线偏微分方程数值解法(附图及matlab程序)

一维抛物线偏微分方程数值解法(4) 上一篇参看一维抛物线偏微分方程数值解法(3)(附图及matlab程序) 解一维抛物线型方程(理论书籍可以参看孙志忠:偏微分方程数值解法) Ut-Uxx=0, 00) U(x,0)=e^x, 0<=x<=1, U(0,t)=e^t,U(1,t)=e^(1+t), 0

一维抛物线偏微分方程数值解法(3)(附图及matlab程序)

一维抛物线偏微分方程数值解法(3) 上一篇参看一维抛物线偏微分方程数值解法(2)(附图及matlab程序) 解一维抛物线型方程(理论书籍可以参看孙志忠:偏微分方程数值解法) Ut-Uxx=0, 00) U(x,0)=e^x, 0<=x<=1, U(0,t)=e^t,U(1,t)=e^(1+t), 0

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软件上的实现张体强1026222 廖荣发1026226 [摘要] 数学物理方程的数值解在实际生活中越来越使用,首先基于偏微分数值解的思想上,通过matlab软件的功能,研究其数学物理方程的数值解,并通过对精确解和数值解进行对比,追究其数值解的可行性,在此,给出相关例子和程序代码,利于以后的再次研究和直接使用。 [关键字] 偏微分方程数值解matlab 数学物理方程的可视化 一:研究意义 在我们解数学物理方程,理论上求数学物理方程的定解有着多种解法,但是有许多定解问题却不能严格求解,只能用数值方法求出满足实际需要的近似解。而且实际问题往往很复杂,这时即便要解出精确解就很困难,有时甚至不可能,另一方面,在建立数学模型时,我们已作了很多近似,所以求出的精确解也知识推导出的数学问题的精确解,并非真正实际问题的精确解。因此,我们有必要研究近似解法,只要使所求得的近似解与精确解之间的误差在规定的范围内,则仍能满足实际的需要,有限差分法和有限元法是两种最常用的

求解数学物理方程的数值解法,而MATLAB 在这一方面具有超强的数学功能,可以用来求其解。 二:数值解法思想和步骤 2.1:网格剖分 为了用差分方法求解上述问题,将求解区域 {}(,)|01,01x t x t Ω=≤≤≤≤作剖分。将空间区间[0,1]作m 等分,将时 间[0,1]区间作n 等分,并记 1/,1/,,0,,0j k h m n x jh j m t k k n ττ===≤≤=≤≤。分别称h 和τ 为空间和 时间步长。用两簇平行直线,0,,0j k x x j m t t k n =≤≤=≤≤将Ω分割成矩形网格。 2.2:差分格式的建立 0u u t x ??-=??………………………………(1) 设G 是,x t 平面任一有界域,据Green 公式(参考数学物理方程第五章): ( )()G u u dxdt udt udx t x Γ??-=--??? ? 其中G Γ=?。于是可将(1)式写成积分守恒形式: ()0udt udx Γ --=? (2) 我们先从(2)式出发构造熟知的Lax 格式设网格如下图所示

Matlab关于数值计算的实现

Matlab关于数值计算的实现 摘要:数值计算(numerical computation computation),主要研究更好的利用计算机更好的进行数值计算,解决各种数学问题。数值分析包括离散傅里叶变换,考虑截断误差,计算误差,函数的敛散性与稳定性等。在数学方面,数值计算的主要研究数值微分与积分,数据的处理与多项式计算,最优化问题,线性方程与非线性方程的求解,常微分方程的数值求解等。同时,数值计算在物理,化学,经济等方面也有研究,本文暂且不表。M atlab软件历经二十多年来的发展,已成为风靡世界的数学三大软件(matlb,Mathematica l,Maple)之一,在数学类科技应用软件中在数值计算方面首屈一指。Matlab以矩阵为数据操作的基本单位,使得矩阵运算十分便捷快速,同时Matlab还提供了海量的计算函数,而且使用可靠地算法进行计算,能使用户在繁复的数学运算中解脱,Matlab还具有方便且完善的图形处理功能,方便绘制二维和三维图形并修饰。

目录 1.数值分析(离散傅里叶变换,考虑截断误差,计算误差,函数 的敛散性与稳定性) 2.数值计算(数值微分与积分,数据的处理与多项式计算, 最优化问题,线性方程与非线性方程的求解,常微分方程的数值求解) 3.图形处理功能(方便绘制二维和三维图形并修饰) 4.总结

1.数据统计与分析 Matlab 可以进行求矩阵的最大最小元素,平均值与中值,关于矩阵元素的求和与求积,累加和与累乘积,标准方程,相关系数,元素排序。现在以求标准方差举例说明Matlab 的实现。 在Matlab 中,实现标准方差计算的函数为std 。对于向量(Y ),std (Y )实现返回一个标准方差,而对于矩阵(A ),std (A )返回一个行向量,该行向量的每个元素对应着矩阵A 各行或各列的标准方差。一般调用std 函数的格式为std (A ,flag ,dim ) Dim 取1或者2分别对应求各列或各行的标准方差,flag 取1时,按照标准方差的计算公式 ∑-=-=N i x x S i N 1 2 1)(11来计算。若flag 取2,则用公式 ∑-==N i x x S i N 1 2 2) (1 进行计算。默认的flag 取值为0,dim 取值为1。课本page143 2. 离散傅里叶变换 离散傅里叶变换广泛应用于信号的分析,光谱和声谱分析、全息技术等各个领域。但直接计算dft 的运算量与变化的长度N 的平方成正比,当N 较大时,计算量太大。随着计算机技术的迅速发展,在计算机上进行离散傅里叶变换计算成为可能。特别是快速傅里叶变换算法的出现,为傅里叶变换的应创造了条件。 (1):傅里叶变换算法的简述。 傅立叶变换是一种分析信号的方法,它可分析信号的成分,也可用这些成分合成信号。许多波形可作为信号的成分,比如正弦波、方波、锯齿波等,傅立叶变换用正弦波作为信号的成分. f(t)是t 的周期函数,如果t 满足狄里赫莱条件:在一个以2T 为周期内f(X)连续或只有有限个第一类间断点,附f (x )单调或可划分成有限个单调区间,则F (x )以2T 为周期的傅里叶级数收敛,和函数S (x )也是以2T 为周期的周期函数,且在这些间断点上,函数是有限值;在一个周期内具有有限个极值点;绝对可积。则有下图①式成立。称为积分运算f(t)的傅立叶变换, ②式的积分运算叫做F(ω)的傅立叶逆变换。F(ω)叫做f(t)的像函数,f(t)叫做 F(ω)的像原函数。F(ω)是f(t)的像。f(t)是F(ω)原像。 ①傅立叶变换 ②傅立叶逆变换

数学应用软件作业6-用Matlab求解微分方程(组)的解析解和数值解

数学应用软件作业6-用Matlab 求解微分方程(组)的解析解和数值解

注:上机作业文件夹以自己的班级姓名学号命名,文件夹包括如下上机报告和Matlab程序。 上机报告模板如下: 佛山科学技术学院 上机报告 课程名称数学应用软件 上机项目用Matlab求解微分方程(组)的解析解和数值解 专业班级姓名学号 一. 上机目的 1.了解求微分方程(组)的解的知识。 2.学习Matlab中求微分方程的各种解的函数,如 dsolve命令、ode45函数等等,其中注意把方程化为新的方程的形式。 3.掌握用matlab编写程序解决求解微分方程的问 题。 二. 上机内容 1、求高阶线性齐次方程:y’’’-y’’-3y’+2y=0。 2、求常微分方程组

2 210cos,2 24,0 t t t dx dy x t x dt dt dx dy y e y dt dt = - = ? +-== ?? ? ?++== ?? 3、求解 分别用函数ode45和ode15s计算求解,分别画出图形,图形分别标注标题。 4、求解微分方程 ,1 )0( ,1 '= + + - =y t y y 先求解析解,在[0,1]上作图; 再用ode45求数值解(作图的图形用“o”表示),在同一副图中作图进行比较,用不同的颜色表示。 三. 上机方法与步骤 给出相应的问题分析及求解方法,并写出Matlab 程序,并有上机程序显示截图。 题1:直接用命令dsolve求解出微分方程的通解。 Matlab程序:

dsolve('D3y-D2y-3*Dy+2*y','x') 题2:将微分方程组改写为 5cos2exp(2) 5cos2exp(2) (0)2,(0)0 dx t t x y xt dy t t x y dt x y ? =+--- ? ? ? =-+-+- ? ? == ? ? ? , 再用命令dsolve求解微分方程的通解。 Matlab程序: 建立timu2.m如下: [x,y]=dsolve('Dx=5*cos(t)+2*exp(-2*t)-x-y','Dy=-5*cos(t)+2*exp(-2*t)+x-y ','x(0)=2,y(0)=0','t') x=simple(x) y=simple(y)

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编写S函数求解微分方程

利用matlab编写S函数求解微分方程自动化专业综合设计报告 自动化专业综合设计报告

函数求解微S编写设计题目:利用 matlab 分方程 自动化系统仿真实验室所在 实验室: 郭卫平 指导教师: 律迪迪学生姓名 200990519114 班级文自0921 学号 成绩评定: 自动化专业综合设计报告

一、设计目的 了解使用simulink的扩展工具——S-函数,s函数可以利用matlab的丰富资源,而不仅仅局限于simulink提供的模块,而用c或c++等语言写的s函数还可以实现对硬件端口的操作,还可以操作windows API 等的,它的魅力在于完美结合了simulink 框图简洁明快的特点和编程灵活方便的优点,提供了增强和扩展sinulink能力的强大机制,同时也是使用RTW实现实时仿真的关键。 二、设计要求 求解解微分方程 y'=y-2x/y 自动化专业综合设计报告 y(0)=1 要求利用matlab编写S函数求解 三、设计内容(可加附页) 【步骤1】获取状态空间表达式。

在matlab中输入 dsolve(‘Dy=y-2*x/y','y(0)=1', 'x') 得到 y=(2*x+1).^(1/2); 【步骤2】建立s函数的m文件。 利用21·用S函数模板文件。 以下是修改之后的模板文件sfuntmpl.m 的内容。 function [sys,x0,str,ts] = sfuntmpl(t,x,u,flag) %SFUNTMPL S-function M-file General template define you can With % M-file S-functions, you own ordinary differential system equations (ODEs), discrete % equations, and/or just about any type of algorithm to be used within a %

一维抛物线型方程数值解法(1)(附图及matlab程序)

一维抛物线偏微分方程数值解法(1) 解一维抛物线型方程(理论书籍可以参看孙志忠:偏微分方程数值解法) Ut-Uxx=0, 00) U(x,0)=e^x, 0<=x<=1, U(0,t)=e^t,U(1,t)=e^(1+t), 00) %不用解线性方程组,由下一层(时间层)的值就直接得到上一层的值 %m,n为x,t方向的网格数,例如(2-0)/0.01=200; %e为误差,p为精确解 u=zeros(n+1,m+1); x=0+(0:m)*h1; t=0+(0:n)*h2; for(i=1:n+1) u(i,1)=exp(t(i)); u(i,m+1)=exp(1+t(i)); end for(i=1:m+1) u(1,i)=exp(x(i)); end for(i=1:n+1) for(j=1:m+1) f(i,j)=0; end end r=h2/(h1*h1); %此处r=a*h2/(h1*h1);a=1 要求r<=1/2差分格式才稳定for(i=1:n) for(j=2:m) u(i+1,j)=(1-2*r)*u(i,j)+r*(u(i,j-1)+u(i,j+1))+h2*f(i,j); end end for(i=1:n+1) for(j=1:m+1) p(i,j)=exp(x(j)+t(i)); e(i,j)=abs(u(i,j)-p(i,j)); end end

数学模型之微分方程及其MATLAB求解

数学模型之微分方程及其MATLAB求解 ---卫星轨迹等经典例题求解分析1. 考虑初值问题画图 y'''?3y ''?y 'y = 0 y(0) = 0 y '(0) =1 y ' '(0) = ?1 2、 3、 【实验步骤与程序】 1. M -文件建立m函数文件

function y=f(t,x) y=[x(2);x(3);9*x(3)^2+x(1)*x(2)]; 求解微分方程,命令如下: x0=[0;1;-1]; [t,y]=ode45(@mm,[0,2.5],x0); plot(y(:,1),y(:,2)); figure(2); plot3(y(:,1),y(:,2),y(:,3))

2、M -文件建立m函数文件 function dx=appollo(t,x) mu=1/82.45; mustar=1-mu; r1=sqrt((x(1)+mu)^2+x(3)^2); r2=sqrt((x(1)-mustar)^2+x(3)^2); dx=[x(2) 2*x(4)+x(1)-mustar*(x(1)+mu)/r1^3-mu*(x(1)-mustar)/r2^3 x(4) -2*x(2)+x(3)-mustar*x(3)/r1^3-mu*x(3)/r2^3];

求解微分方程,命令如下: x0=[1.2;0;0;-1.04935751]; options=odeset('reltol',1e-8); [t,y]=ode45(@appollo,[0,20],x0,options); plot(y(:,1),y(:,3)) title('Appollo卫星运动轨迹') xlabel('x') ylabel('y')

双曲方程基于matlab的数值解法

双曲型方程基于MATLAB 的数值解法 (数学1201,陈晓云,41262022) 一:一阶双曲型微分方程的初边值问题 0,01,0 1.(,0)cos(),0 1. (0,)(1,)cos(),0 1. u u x t t x u x x x u t u t t t ππ??-=≤≤≤≤??=≤≤=-=≤≤ 精确解为 ()t x cos +π 二:数值解法思想和步骤 2.1:网格剖分 为了用差分方法求解上述问题,将求解区域{}(,)|01,01x t x t Ω=≤≤≤≤作剖分。将空间区间[0,1]作m 等分,将时间[0,1]区间作n 等分,并记 1/,1/,,0,,0j k h m n x jh j m t k k n ττ===≤≤=≤≤。分别称h 和τ为空间和时 间步长。用两簇平行直线,0,,0j k x x j m t t k n =≤≤=≤≤将Ω分割成矩形网格。 2.2:差分格式的建立 0u u t x ??-=?? 2.2.1:Lax-Friedrichs 方法 对时间、空间采用中心差分使得 2h 1 1111)(2 1u u x u u u u u t u k j k j k j k j k j k j -+-++-= +=-= ????τ τ 则由上式得到Lax-Friedrichs 格式 1 11111()202k k k k k j j j j j u u u u u h τ+-+-+-+-+=

截断误差为 ()[]k k k j h j j R u L u Lu =- 1 11111()22k k k k k k k j j j j j j j u u u u u u u h t x τ+-+-+-+-??=+-+?? 23222 3 (),(0,0)26k k j j u u h O h j m k n t x ττ??= -=+≤≤≤≤?? 所以Lax-Friedrichs 格式的截断误差的阶式2()O h τ+ 令/s h τ=:则可得差分格式为 1111 11(),(0,0)222 k k k k k j j j j j s s u u u u u j m k n +--++=-+++≤≤≤≤ 0cos()(0)j j u x j m π=≤≤ 0cos(),cos(),(0)k k k m k u t u t k n ππ==-≤≤ 其传播因子为: ()()()e e G h i h i s h i h i σσσστσ---=-+e e 221, 化简可得: ()()()()()h s G h is h G στσσστ σsin 11,sin cos ,2 2 2--=-= 所以当1s ≤时,()1,≤τσG ,格式稳定。 * 2.2.2:LaxWendroff 方法 用牛顿二次插值公式可以得到LaxWendroff 的差分格式,在此不详细分析,它的截断误差为() h 2 2 +O τ ,是二阶精度;当2s ≤时,()1,≤τσG , 格式稳定。在这里主要用它与上面一阶精度的Lax-Friedrichs 方法进行简单对比。 2.3差分格式的求解

一维偏微分方程的pdepe(matlab)函数 解法

本文根据matlab帮助进行加工,根据matlab帮助上的例子,帮助更好的理解一维偏微分方程的pdepe函数解法,主要加工在于程序的注释上。 Examples Example 1.This example illustrates the straightforward formulation, computation, and plotting of the solution of a single PDE. This equation holds on an interval for times . The PDE satisfies the initial condition and boundary conditions It is convenient to use subfunctions to place all the functions required by pdepe in a single function. function pdex1 m = 0; x = linspace(0,1,20); %linspace(x1,x2,N)linspace是Matlab中的一个指令,用于产生x1,x2之间的N点行矢量。 %其中x1、x2、N分别为起始值、终止值、元素个数。若缺省N,默认点数为100 t = linspace(0,2,5); sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);

% Extract the first solution component as u. u = sol(:,:,1); % A surface plot is often a good way to study a solution. surf(x,t,u) title('Numerical solution computed with 20 mesh points.') xlabel('Distance x') ylabel('Time t') % A solution profile can also be illuminating. figure plot(x,u(end,:)) title('Solution at t = 2') xlabel('Distance x') ylabel('u(x,2)') % -------------------------------------------------------------- function [c,f,s] = pdex1pde(x,t,u,DuDx) c = pi^2; f = DuDx; s = 0; % -------------------------------------------------------------- function u0 = pdex1ic(x) u0 = sin(pi*x); % -------------------------------------------------------------- function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = pi * exp(-t); qr = 1;

相关主题
文本预览
相关文档 最新文档