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类似，但它的

符号计算不是基于Maple上的，而是自己开发的。

Mathematica的基本系统主要是用C语言开发的，因而可以比较容易地移植到各种平台上，Mathematica是一个交互式的计算系统，计算是在用户和Mathematica互相交换、传递信息数据的过程中完成的。Mathematica系统所接受的命令都被称作表达式，系统在接受了一个表达式之后就对它进行处理，然后再把计算结果返回。Mathematica对于输入形式有比较严格的规定，用户必须按照系统规定的数学格式输入，系统才能正确地处理，不过由于3. 0版本引入输入面板，并且可以修改、重组输入面板，因此以前版本输入指令时需要不断切换大小写字符的繁琐方式得到很好的改善。3.0版本可以用各种格式保存文件和剪贴内容，包括RTF、HTML、BMP等格式。

五、四种软件的比较

选用何种数学软件?如果仅仅是要求一般的计算或者是普通用户日常使用，首选的是M athCAD，它在高等数学方面所具有的能力，足够一般客户的要求，而且它的输入界面也特别友好。如果要求计算精度、符号计算和编程方面的话，最好同时使用Maple和Mathemat ica，它们在符号处理方面各具特色，有些Maple不能处理的，Mathematica却能处理，诸如某些积分、求极限等方面，这些都是比较特殊的。如果要求进行矩阵方面或图形方面的处理，则选择MATLAB，它的矩阵计算和图形处理方面则是它的强项，同时利用MATLAB的N oteBook功能，结合Word6.0/7.0的编辑功能，可以很方便地处理科技文章。

1.2 满足Neumann边界条件的Helmholtz方程

源程序：

function [u,x,y] = Helmholtz_Newton(f,g,dbx,bx0,bxf,by0,byf,D,Mx,My,MinErr,M axIter)

%解方程：u_xx + u_yy + g(x,y)u = f(x,y)

% 自变量取值区域D = [x0,xf,y0,yf] = {(x,y) |x0 <= x <= xf, y0 <= y <= y f}

% 边界条件

% u(x0,y) = bx0(y), u(xf,y) = bxf(y)

% u(x,y0) = by0(x), u(x,yf) = byf(x)

% x轴均分为Mx段

% y轴均分为My段

% tol 误差因子

% MaxIter: 最大迭代次数

x0 = D(1); xf = D(2); y0 = D(3); yf = D(4);

dx = (xf - x0)/Mx; x = x0 + [0:Mx]*dx;%构造内点数组

dy = (yf - y0)/My; y = y0 + [0:My]'*dy;

Mx1 = Mx + 1; My1 = My + 1;

for i = 1:Mx

for j = 1:My

F(i,j) = f(x(i),y(j)); G(i,j) = g(x(i),y(j));

end

end

dx2 = dx*dx; dy2 = dy*dy; dxy2 = 2*(dx2 + dy2);

rx = dx2/dxy2; ry = dy2/dxy2; rxy = rx*dy2;

%边界条件

for m = 1:My1

u([1 Mx1],m)=[bx0(y(m)) bxf(y(m))]; %左右边界

end

for n = 1:Mx1

u(n,[1 My1]) = [by0(x(n)); byf(x(n))];%上下边界

end

%边界平均值作迭代初值

sum_of_bv = sum(sum([u([1 Mx1],2:My) u(2:Mx,[1 My1])']));

u(2:Mx,2:My) = sum_of_bv/(2*(Mx + My - 2));

for itr = 1:MaxIter

for i = 2:(Mx1-1)

u(i,1)=2*ry*u(i,2)+rx*(u(i+1,1)+u(i-1,1))+rxy*(G(i,1)*u(i,1)-F(i,1)-2*dbx (x(i),y(1))/dx);

end

for j=2:(My1-1)

u(1,j)=ry*(u(1,j+1)+u(1,j-1))+2*rx*u(2,j)+rxy*(G(1,j)*u(1,j)-F(1,j)-2*dbx (x(1),y(j))/dy);

end

for i= 2:Mx

for j = 2:My

u(i,j) = ry*(u(i+1,j)+u(i-1,j)) + rx*(u(i,j-1)+u(i,j+1))+ rxy*(G(i,j)*u (i,j)- F(i,j)); %迭代公式

end

end

if itr > 1 & max(max(abs(u - u0))) < MinErr%循环结束条件

break;

end

u0 = u;

end

u=u';

例1.1迭代法求解满足Neumann型边界条件的Helmholtz方程应用实例。求以下满足Ne umann型边界条件的的数值解：

自变量取值：

边界:

解：可知，

在MATLAB中编写脚本文件：

f = inline('x^2+y^2','x','y');

g = inline('sqrt(x)','x','y');

x0 = 0; xf = 4; y0 = 0; yf = 4;%自变量取值范围

Mx = 50;My = 30;%等分段数

dbx=inline('x^2+y^2','x','y');

bx0 = inline('y^2','y'); %边界条件

bxf = inline('4^2*cos(y)','y');

by0 = inline('x^2','x');

byf = inline('4^2*cos(x)','x');

D = [x0 xf y0 yf]; MaxIter = 100; MinErr = 1e-4;

[U,x,y] = Helmholtz_Newton(f,g,dbx,bx0,bxf,by0,byf,D,Mx,My,MinErr,MaxIter); clf, mesh(U)

xlabel('x')

ylabel('y')

zlabel('U')

2 抛物形偏微分方程

2.1显式前向欧拉法

源程序：

function [u,x,t] = EF_Euler(A,xf,T,it0,bx0,bxf,M,N)

%解方程A u_xx = u_t ，0 <= x <= xf, 0 <= t <= T

%初值: u(x,0) = it0(x)

% 边界条件: u(0,t) = bx0(t), u(xf,t) = bxf(t)

% M ：x 轴的等分段数

% N ：t 轴的等分段数

dx = xf/M; x = [0:M]*dx;

dt = T/N; t = [0:N]'*dt;

for i= 1:M + 1

u(i,1) = it0(x(i));

end

for j = 1:N + 1

u([1 M + 1],j) = [bx0(t(j)); bxf(t(j))];

end

r = A*dt/dx/dx, r1 = 1 - 2*r;

if(r>0.5)

disp('r>0.5,unstability');

end

for j = 1:N

for i = 2:M

u(i,j+1) = r*(u(i + 1,j) + u(i-1,j)) + r1*u(i,j); %(9.2.3)

end

end

u=u';

例2.1显式前向欧拉法求解一维抛物性方程应用实例。求满足以下条件的热传导数值解：

自变量取值：

边界:

解：在MATLAB中编写脚本文件：

A = 0.5; %方程系数

it0 = inline('sin(pi*x)','x'); %初始条件

bx0 = inline('0'); bxf = inline('0'); %边界条件

xf = 2; M = 80; T = 0.1; N = 100;

[u1,x,t] = EF_Euler(A,xf,T,it0,bx0,bxf,M,N);

figure(1),clf,mesh(u1)

xlabel('x')

ylabel('t')

zlabel('U')

title('r>0.5')

M=50;

[u1,x,t] = EF_Euler(A,xf,T,it0,bx0,bxf,M,N);

figure(2),clf,mesh(u1)

xlabel('x')

ylabel('t')

zlabel('U')

title('r<0.5')

2.2 隐式后向欧拉法

源程序：

function [u,x,t] = IB_Euler(A,xf,T,it0,bx0,bxf,M,N)

%解方程A1 u_xx = u_t ，0 <= x <= xf, 0 <= t <= T

%初值: u(x,0) = it0(x)

% 边界条件: u(0,t) = bx0(t), u(xf,t) = bxf(t)

% M ：x 轴的等分段数

% N ：t 轴的等分段数

dx = xf/M; x = [0:M]*dx;

dt = T/N; t = [0:N]'*dt;

for i = 1:M + 1

u(i,1) = it0(x(i));

end

for j = 1:N + 1

u([1 M + 1],j) = [bx0(t(j)); bxf(t(j))];

end

r = A*dt/dx/dx;

r2 = 1 + 2*r;

for i = 1:M - 1

P(i,i) = r2; %构造9.2.9的矩阵

if i > 1

P(i - 1,i) = -r; P(i,i - 1) = -r;

end

end

for j = 2:N + 1

b = [r*u(1,j); zeros(M - 3,1); r*u(M + 1,j)] + u(2:M,j - 1); %Eq.(9.2.9)

u(2:M,j) = linsolve(P,b);

end

u=u';

例2.1隐式前向欧拉法求解一维抛物性方程应用实例。求满足以下条件的热传导数值解：

自变量取值：

边界:

解：在MATLAB中编写脚本文件：

A= 0.5; %方程系数

it0 = inline('sin(pi*x)','x'); %初始条件

bx0 = inline('0'); bxf = inline('0'); %边界条件

xf = 2; M = 50; T = 0.1; N = 100;

[u1,x,t] = IB_Euler(A,xf,T,it0,bx0,bxf,M,N);

mesh(u1)

xlabel('x')

ylabel('t')

zlabel('U')

2.3 Grank-Nicholson方法

源程序：

function [u,x,t] = Grank_Nicholson(A,xf,T,it0,bx0,bxf,M,N)

%解方程A u_xx = u_t ，0 <= x <= xf, 0 <= t <= T

%初值: u(x,0) = it0(x)

% 边界条件: u(0,t) = bx0(t), u(xf,t) = bxf(t)

% M ：x 轴的等分段数

% N ：t 轴的等分段数

dx = xf/M; x = [0:M]*dx;

dt = T/N; t = [0:N]'*dt;

for i = 1:M + 1

u(i,1) = it0(x(i));

end

for n = 1:N + 1

u([1 M + 1],n) = [bx0(t(n)); bxf(t(n))];

end

r = A*dt/dx/dx;

r1 = 2*(1 + r); r2 = 2*(1 - r);

for i = 1:M - 1

P(i,i) = r1; %(9.2.17)

Q(i,i)=r2;

if i > 1

P(i - 1,i) = -r; P(i,i - 1) = -r; %(9.2.17)等式左边矩阵

Q(i - 1,i) = r; Q(i,i - 1) = r;%(9.2.17)等式右边矩阵

end

end

for k = 2:N + 1

b=Q*u(2:M,k-1)+[r*(u(1,k)+u(1,k-1));zeros(M-2,1)];

u(2:M,k) = linsolve(P,b); %(9.2.17)

end

u=u';

例2.1Grank-Nicholson方法求解一维抛物性方程应用实例。求满足以下条件的热传导数值解：

自变量取值：

边界:

解：在MATLAB中编写脚本文件：

A= 0.5; %方程系数

it0 = inline('sin(pi*x)','x'); %初始条件

bx0 = inline('0'); bxf = inline('0'); %边界条件

xf = 2; M = 25; T = 0.1; N = 100;

[u1,x,t] = Grank_Nicholson(A,xf,T,it0,bx0,bxf,M,N);

mesh(u1)

xlabel('x')

ylabel('t')

zlabel('U')

3双曲线偏微分方程

3.1显式中心差分法

源程序：

function [u,x,t] = ECD_Wave(A,xf,T,it0,i1t0,bx0,bxf,M,N)

%解方程a u_xx = u_tt for 0<=x<=xf, 0<=t<=T

% 初始条件: u(x,0) = it0(x), u_t(x,0) = i1t0(x)

% 边界条件: u(0,t)= bx0(t), u(xf,t) = bxf(t)

% M ：沿x轴的等分段数

% N ：沿y轴的等分段数

dx = xf/M; x = [0:M]'*dx;

dt = T/N; t = [0:N]*dt;

for i = 1:M + 1

u(i,1) = it0(x(i));

end

for k = 1:N + 1

u([1 M + 1],k) = [bx0(t(k)); bxf(t(k))];

end

r = A*(dt/dx)^ 2; r1 = r/2; r2 = 2*(1 - r);

u(2:M,2) = r1*u(1:M - 1,1) + (1 - r)*u(2:M,1) + r1*u(3:M + 1,1) + dt*i1t0(x (2:M)); %(11.3.4)

for k = 3:N + 1

u(2:M,k) = r*u(1:M - 1,k - 1) + r2*u(2:M,k-1) + r*u(3:M + 1,k - 1)- u(2: M,k - 2); %(11.3.3)

end

u=u';

例3.1显式中心差分法求解一维波动方程应用实例。求满足以下条件的波动方程数值解：

自变量取值：

边界:

解：在MATLAB中编写脚本文件：

A= 1;

it0 = inline('x-x^2','x'); i1t0 = inline('0');

bx0t = inline('0'); bxft = inline('0');

xf =1; M = 10; T = 2; N = 50;

[u,x,t] = ECD_Wave(A,xf,T,it0,i1t0,bx0t,bxft,M,N);

mesh(u)

xlabel('x')

ylabel('t')

zlabel('U')

（1）通过函数fem_coef为每一个节点和子函数构造基函数。

function [U,c] = fem_coef(f,g,p,c,N,S,N_i)

%p(i,s,1:3): 基函数ftn phi_i系数

%c = [ .1 1 . 0 0 .] 边界节点取1，内节点取0

%N(n,1:2) : 第n个节点的x和y坐标

%S(s,1:3) : 第s个子区域的节点#s

%N_i : 内节点个数

%U(s,1:3) : 每个区域的p1 + p2(s)x + p3(s)y 的坐标

N_n = size(N,1); % 总共节点数

N_s = size(S,1); % 总共子区域数

d=zeros(N_i,1);

N_b = N_n-N_i;

for i = N_b+1:N_n

for n = 1:N_n

for s = 1:N_s

xy = (N(S(s,1),:) + N(S(s,2),:) + N(S(s,3),:))/3; %重心

p_vctr = [p([i n],s,1) p([i n],s,2) p([i n],s,3)];

tmpg(s) = sum(p(i,s,2:3).*p(n,s,2:3))-g(xy(1),xy(2))*p_vctr(1,:)*[1 x y]'*p_vctr(2,:)*[1 xy]';

dS(s) = det([N(S(s,1),:) 1; N(S(s,2),:) 1;N(S(s,3),:) 1])/2;

%子区域

if n == 1, tmpf(s) = -f(xy(1),xy(2))*p_vctr(1,:)*[1 xy]'; end end

A12(i - N_b,n) = tmpg*abs(dS)';

end

d(i-N_b) = tmpf*abs(dS)';

end

d = d - A12(1:N_i,1:N_b)*c(1:N_b)';

c(N_b + 1:N_n) = A12(1:N_i,N_b+1:N_n)\d;

for s = 1:N_s

for j = 1:3, U(s,j) = c*p(:,s,j); end

end

（2）得到基函数图像

function Show_Basis()

%N = [-1 1;1 1;1 -1;-1 -1;0.2 0.5]; %节点集合

N = [1 0;0 1;1 2;2 1;1.2 1.5]; %节点集合

N_n = size(N,1); % 总节点数

S = [1 2 5;2 3 5;3 4 5;1 4 5]; %区域集合

N_s = size(S,1); % 总区域数

figure(1), clf

for s = 1:N_s

nodes = [S(s,:) S(s,1)];

for i = 1:3

plot([N(nodes(i),1) N(nodes(i + 1),1)],[N(nodes(i),2) N(nodes(i+1),2)]), hold on

end

end

text(0.8,0.6,'S1');text(0.8,1.6,'S2');text(1.4,1.5,'S3'),text(1.4,0.6,'S4');

%基函数

p = fem_basis_ftn(N,S);

%x0 = -1; xf = 1; y0 = -1; yf = 1; %graphic region

x0 = 0; xf = 2; y0 = 0; yf = 2; %graphic region

figure(2), clf

Mx = 50; My = 50;

c = [0 1 2 3 0]; %节点的值

dx = (xf - x0)/Mx; dy = (yf - y0)/My;

xi = x0 + [0:Mx]*dx; yi = y0 + [0:My]*dy;

i_ns = [1 2 3 4 5]; %节点标号

for itr = 1:5

i_n = i_ns(itr);

if itr == 1

for i = 1:length(xi)

for j = 1:length(yi)

Z(j,i) = 0;

for s = 1:N_s

if inpolygon(xi(i),yi(j), N(S(s,:),1),N(S(s,:),2)) > 0

Z(j,i) = p(i_n,s,1) + p(i_n,s,2)*xi(i) + p(i_n,s,3)*yi (j);

break;

end

end

end

end

subplot(321), mesh(xi,yi,Z);title(itr) %节点1的基函数

else

c1 = zeros(size(c)); c1(i_n) = 1;

subplot(320 + itr)

trimesh(S,N(:,1),N(:,2),c1) %节点2-5的基函数

title(itr);

end

end

c = [0 1 2 3 0];

subplot(326)

trimesh(S,N(:,1),N(:,2),c);title('基函数的线型组合')