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

1、Maple、MATLAB、MathCAD和Mathematica 2008-04-25 11:43 | (分类:默认分类) 一、Maple V 系统Maple V是由Waterloo大学开发的数学系统软件，它不但具有精确的数值处理功能，而且具有无以伦比的符号计算功能。Maple V的符号计算能力还是MathCAD和MATLAB等软件的符号处理的核心。Maple提供了2000余种数学函数，涉及范围包括：普通数学、高等数学、线性代数、数论、离散数学、图形学。它还提供了一套内置的编程语言，用户可以开发自己的应用程序，而且Maple自身的2000多种函数，基本上是用此语言开发的。 Maple采用字符行

2、输入方式，输入时需要按照规定的格式输入，虽然与一般常见的数学格式不同，但灵活方便，也很容易理解。输出则可以选择字符方式和图形方式，产生的图形结果可以很方便地剪贴到Windows应用程序内。 二、MATLAB 系统 MATLAB原是矩阵实验室（Matrix Laboratory）在70年代用来提供Linpack和Eispack软件包的接口程序，采用C语言编写。从80年代出现3.0的DOS版本，逐渐成为科技计算、视图交互系统和程序语言。MATLAB可以运行在十几个操作平台上，比较常见的有基于Windows 9X/NT、OS/2、Macintosh、Sun、Unix、Linux等平台的系统。 MAT

3、LAB程序主要由主程序和各种工具包组成，其中主程序包含数百个内部核心函数，工具包则包括复杂系统仿真、信号处理工具包、系统识别工具包、优化工具包、神经网络工具包、控制系统工具包、分析和综合工具包、样条工具包、符号数学工具包、图像处理工具包、统计工具包等。而且5.x版本还包含一套几十个的PDF文件，从MATLAB的使用入门到其他专题应用均有详细的介绍。 MATLAB是数值计算的先锋，它以矩阵作为基本数据单位，在应用线性代数、数理统计、自动控制、数字信号处理、动态系统仿真方面已经成为首选工具，同时也是科研工作人员和大学生、研究生进行科学研究的得力工具。MATLAB在输入方面也很方便，可以使用内部的E

4、ditor或者其他任何字符处理器，同时它还可以与Word6.0/7.0结合在一起，在Word的页面里直接调用MATLAB的大部分功能，使Word具有特殊的计算能力。 三、MathCAD 系统 MathCAD是美国Mathsoft公司推出的一个交互式的数学系统软件。从早期的DOS下的1.0和Windows下的4.0版本，到今日的8.0版本，功能也从简单的数值计算，直至引用Maple强大的符号计算能力，使得它发生了一个质的飞跃。 MathCAD是集文本编辑、数学计算、程序编辑和仿真于一体的软件。MathCAD7.0 Professional（专业版）运行在Win9X/NT下，它的主要特点是输入格式

5、与人们习惯的数学书写格式很近似，采用WYSWYG（所见所得）界面，特别适合一般无须进行复杂编程或要求比较特殊的计算。MathCAD 7.0 Professional 还带有一个程序编辑器，对于一般比较短小，或者要求计算速度比较低时，采用它也是可以的。这个程序编辑器的优点是语法特别简单。MathCAD可以看作是一个功能强大的计算器，没有很复杂的规则；同时它也可以和Word、Lotus、WPS2000等字处理软件很好地配合使用，可以把它当作一个出色的全屏幕数学公式编辑器。 四、Mathematica 系统 Mathematica是由美国物理学家Stephen Wolfram领导的Wolfram R

6、esearch开发的数学系统软件。它拥有强大的数值计算和符号计算能力，在这一方面与Maple类似，但它的符号计算不是基于Maple上的，而是自己开发的。 Mathematica的基本系统主要是用C语言开发的，因而可以比较容易地移植到各种平台上，Mathematica是一个交互式的计算系统，计算是在用户和Mathematica互相交换、传递信息数据的过程中完成的。Mathematica系统所接受的命令都被称作表达式，系统在接受了一个表达式之后就对它进行处理，然后再把计算结果返回。Mathematica对于输入形式有比较严格的规定，用户必须按照系统规定的数学格式输入，系统才能正确地处理，不过由于3

7、.0版本引入输入面板，并且可以修改、重组输入面板，因此以前版本输入指令时需要不断切换大小写字符的繁琐方式得到很好的改善。3.0版本可以用各种格式保存文件和剪贴内容，包括RTF、HTML、BMP等格式。 五、四种软件的比较 选用何种数学软件?如果仅仅是要求一般的计算或者是普通用户日常使用，首选的是MathCAD，它在高等数学方面所具有的能力，足够一般客户的要求，而且它的输入界面也特别友好。如果要求计算精度、符号计算和编程方面的话，最好同时使用Maple和Mathematica，它们在符号处理方面各具特色，有些Maple不能处理的，Mathematica却能处理，诸如某些积分、求极限等方面，这些都

8、是比较特殊的。如果要求进行矩阵方面或图形方面的处理，则选择MATLAB，它的矩阵计算和图形处理方面则是它的强项，同时利用MATLAB的NoteBook功能，结合Word6.0/7.0的编辑功能，可以很方便地处理科技文章。MATLAB精通科学计算之十偏微分方程求解（二） 2008-09-14 09:29 | (分类:默认分类) 1.2 满足Neumann边界条件的Helmholtz方程源程序：function u,x,y = Helmholtz_Newton(f,g,dbx,bx0,bxf,by0,byf,D,Mx,My,MinErr,MaxIter)%解方程： u_xx + u_yy + g(

9、x,y)u = f(x,y)% 自变量取值区域 D = x0,xf,y0,yf = (x,y) |x0 = x = xf, y0 = y 1 & max(max(abs(u - u0) MinErr%循环结束条件 break; end u0 = u;endu=u;例1.1 迭代法求解满足Neumann型边界条件的Helmholtz方程应用实例。求以下满足Neumann型边界条件的的数值解：自变量取值：边界:解：可知，在MATLAB中编写脚本文件：f = inline(x2+y2,x,y); g = inline(sqrt(x),x,y);x0 = 0; xf = 4; y0 = 0; yf =

10、 4;%自变量取值范围Mx = 50;My = 30;%等分段数dbx=inline(x2+y2,x,y);bx0 = inline(y2,y); %边界条件bxf = inline(42*cos(y),y); by0 = inline(x2,x); byf = inline(42*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(

11、y)zlabel(U)MATLAB精通科学计算之十偏微分方程求解（三） 2008-09-14 09:30 | (分类:默认分类) 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 0.5) disp(r0.5,unstability);endfor 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) endendu=u;例

12、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(r0.5)M=50;u1,x,t =

13、 EF_Euler(A,xf,T,it0,bx0,bxf,M,N);figure(2),clf,mesh(u1)xlabel(x)ylabel(t)zlabel(U)title(r0.5)MATLAB精通科学计算之十偏微分方程求解（四） 2008-09-14 09:30 | (分类:默认分类) 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 1 P(i - 1,i) = -r; P(i,i - 1) = -r; endendfor j =

14、 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);endu=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,

15、x,t = IB_Euler(A,xf,T,it0,bx0,bxf,M,N); mesh(u1)xlabel(x)ylabel(t)zlabel(U)MATLAB精通科学计算之十偏微分方程求解（五） 2008-09-14 09:31 | (分类:默认分类) 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 1 P(i - 1,i) = -r; P(i,i - 1) = -r; %(9.2.17)等式左边矩阵

16、Q(i - 1,i) = r; Q(i,i - 1) = r;%(9.2.17)等式右边矩阵 endendfor 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)endu=u;例2.1 Grank-Nicholson方法求解一维抛物性方程应用实例。求满足以下条件的热传导数值解： 自变量取值：边界:解：在MATLAB中编写脚本文件：A= 0.5; %方程系数it0 = inline(sin(pi*x),x); %初始条件bx0 = inline(0); b

17、xf = 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)MATLAB精通科学计算之十偏微分方程求解（七） 2008-09-14 09:32 | (分类:默认分类) 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 0 Z(j,i) = p(i_n,s,1) + p(i_n,s,2)*xi(i) + p(