malab软件在数值模拟中的应用

malab软件在数值模拟中的应用

1b（magneter）是一个功能强、效率高、科学计算和应用的嵌入式软件包。它是美国Mathworks公司推出的一种高性能的科技计算软件,它集数值计算、符号运算、图形图像处理、编程等多种功能于一体,应用非常广泛。Matlab对于学习者的编程语言基础要求不高,库函数和编程语句丰富多样且简单易学,在数据可视化上也有独特的优势。学习者不需要投入太多的时间在学习编程语言知识上,可以直接利用软件提供的丰富函数,编写较简单的程序即可解决许多普通物理学的问题。一般而言,对从实际问题中建立的方程或公式也很难用手工计算出满意的结果,这就需要有一种方便快捷的方法来使以上问题形象化。基于此,本文探讨Matlab软件在数值模拟方面的应用,为利用Matlab进行计算机模拟提供一些可供借鉴的范例。1以文字生成的方式生成文本Matlab输入命令的方式有两种:一种是在命令窗口中直接输入简单的语句,这种方式适用于命令比较简单且处理问题没有普遍应用性的场合;另一种工作方式为M文件的编程工作方式。M文件是一个简单的文本文件,语法比一般的高级语言都要简单,交互性强;而且可以像一般文本文件那样在任何文本编辑器中进行编辑、存储、修改和读取。当使用绘图语句时,Matlab自动打开一个图形窗口。如果已经有图形窗口存在,作图命令便会使用已存在的图形窗口。如果使用命令figure,就会打开一个新的图形窗口。在已有图形上继续作图的指令是holdon;取消这种功能的指令是holdoff。1.1二维构图的对象1pool属性二维图形绘图命令中最基本的指令就是plot。如果输入两个矢量x、y,则plot(x,y)产生的是y相对于x的图形。2fpot句子当已知函数关系的时候,采用fplot语句可以更加快速、精确地绘制出指定区间上的函数图像。1.23d曲线设计1求取点坐标其中x、y、z是同维数的矢量或矩阵。每组x、y、z构成一个点的坐标,各点依次相连,形成一条曲线。如果是矩阵,则它们相应的列构成一条三维曲线的数据点坐标,所以用矩阵可以同时画多条空间曲线。2生成数据网格的命令Matlab中,曲面是用xy平面的各个格点上的z坐标来定义,相邻点用直线连接。因此,为了描绘二元函数形成的空间曲面,首先要建立平面的数据网格。即根据x、y的定义域[x1,x2]、[y1,y2]生成平面上的坐标格点,每个格点对应一组[x,y]。生成数据网格的命令是meshgrid,格式为[X,Y]=meshgrid(x,y),使用meshgrid命令,将两个独立的矢量x、y组成数据网格。mesh是三维网格作图命令,mesh(x,y,z)画出了每一个格点(x,y)上对应的z值。surfl(X,Y,Z)是画带光照模式的三维曲面图。该命令显示一个带阴影的曲面,融合了周围的、散射的和镜面反射的光照模式。要获得较平滑的颜色过度,应该使用有线性强度变化的色图,如gray、copper、bone、pink等。1.3阶微分方程的求解Matlab解常微分方程组的能力很强且方便,对于在普通物理学中遇到的大多数动力学方程都可用命令ode45求解。Matlab只能解一阶的常微分方程组,高阶的常微分方程需要转化成一阶方程组才能求解。对于二阶常微分方程F(¨x‚˙x‚x,t)=0F(x¨‚x˙‚x,t)=0,首先需要化成显式形式x=f(˙x‚x,t)x=f(x˙‚x,t),然后令y(1)=x‚y(2)=˙xy(1)=x‚y(2)=x˙,则二阶常微分方程化为两个一阶常微分方程组成的方程组dy(1)dt=y(2)dy(1)dt=y(2)和dy(2)dt=f(y(1),y(2)‚t)dy(2)dt=f(y(1),y(2)‚t),从而使问题得到解决。2bp是模拟点电压的电势和车间2.1[x,y,z]的电势及数字图件物理情景是Oxy平面上,在x=0、y=1处有一正电荷,x=0、y=-1处有一负电荷,求出此时的电势分布情况。由电势的表达式U=q4πεrU=q4πεr以及电势叠加原理可以计算出2个点电荷电场的电势分布,其中r=√(x-x0)2+(y-y0)2r=(x−x0)2+(y−y0)2−−−−−−−−−−−−−−−−−√。启动M文件编辑器,创建一个新的M文件,在编辑窗口中键入下列语句:[x,y]=meshgrid(-2∶0.1∶2,-2∶0.1∶2);%建立数据网格z=1./sqrt(x.^2+(y-1).^2+0.01)-1./sqrt(x.^2+(y+1).^2+0.01);%电势的表达式surfl(x,y,z);%三维曲面绘图shadinginterp;%平滑三维曲面colormap(gray);%设定颜色为灰色选定一系列x和y后,组成平面上的网格点,再计算对应每一点上的z值。-2∶0.1∶2,-2∶0.1∶2分别是选取横坐标与纵坐标的一系列数值,meshgrid是生成数据网格的命令。z=1./sqrt(x.^2+(y-1).^2+0.01)-1./sqrt(x.^2+(y+1).^2+0.01)是场点(x,y)的电势。当场点在电荷处时,会出现分母为零的情况,因此,在r里加了一个小量0.01,这样既可以完成计算,又不会对结果的正确性造成太大影响。另外需要注意的是表达式中的“./”、“.^”是对数组运算的算符,含义与数值运算中的“/”、“^”相同,不同之处是后者只对单个数值变量进行运算,而前者对整个数组变量中的所有元素同时进行运算。surfl(x,y,z)函数绘出的三维图形的表面是画带光照模式的三维曲面图,图形输出如图1(a)所示。可以将其换为mesh(x,y,z),其图形输出见图1(b)。2.2电场线的绘制首先建立电场线的微分方程(二维情况)。因为电场中任一点的电场方向都沿该点电场线的切线方向,所以满足dydx=EyEx.dydx=EyEx.引入参变量t得到dxEx=dyEy=t.dxEx=dyEy=t.设二点电荷位于(0,-1)和(0,1),二点电荷电量为q1和q2(均等于1),由库伦定律和电场的叠加原理,得出微分方程˙x=dxdt=Ex=q1x[x2+(y+1)2]3/2+q2x[x2+(y-1)2]3/2‚˙y=dydt=Ey=q1(y+1)[x2+(y+1)2]3/2+q2(y-1)[x2+(y-1)2]3/2.x˙=dxdt=Ex=q1x[x2+(y+1)2]3/2+q2x[x2+(y−1)2]3/2‚y˙=dydt=Ey=q1(y+1)[x2+(y+1)2]3/2+q2(y−1)[x2+(y−1)2]3/2.解此方程就可以绘制出电场线。下面是微分方程的函数文件:functionydot=electricfun(t,y,flag,p1,p2);%p1,p2是参量,表示电量ydot=[p1*y(1)/(sqrt(y(1).^2+(y(2)+1).^2).^3)+p2*y(1)/(sqrt(y(1).^2+(y(2)-1).^2).^3);%dx/dt=Exp1*(y(2)+1)/(sqrt(y(1).^2+(y(2)+1).^2).^3)+p2*(y(2)-1)/(sqrt(y(1).^2+(y(2)-1).^2).^3));%dy/dt=Ey编写好函数文件后,命名为electricfun.m,存在当前路径下,然后开始编写解微分方程的主程序electric.m。p1=1;p2=1;%点电荷所带电量axis([-2,2,-2,2]);%设定坐标轴范围-2≤x≤2,-2≤y≤2holdon;%图形控制,不可擦除模式plot(0,1,′*r′);plot(0,-1,′*r′);%绘制两源电荷a=(pi/24):pi/12:(2*pi-pi/24);%圆周上电场线起点所对应的角度b=0.1*cos(a);c=0.1*sin(a);%电场线起点所对应的相对坐标c1=-1+c;c2=1+c;%把起点圆周的圆心放置在源电荷处b0=[bb];c0=[c1c2];%初始条件,所有电场线的起点横、纵坐标构成了矢量b0和c0fori=1:48;%循环求解48次微分方程[t,y]=ode45(′electricfun′,[0:0.05:40],[b0(i),c0(i)],,p1,p2);%调用ode45求解,对应一个初条件(起点),求解出一条电场线plot(y(:,1),y(:,2),′b′);%绘制出此条电场线end;%结束循环,共绘制出48条电场线在程序中通过从圆周上取24个不同的角度(从π/24到2π-π/24,每隔π/12取1个角度),然后算出每个角度上起点的横、纵坐标值;[b1,c]和[b2,c]分别是以2个源点电荷为圆心