应用MATLAB求解经典物理若干典型问题 毕业论文.docx

本科生毕业论文应用matlab求解经典物理若干典型问题the application of matlab in solving some classical physics questions 专业班级：物理0701学生姓名：*指导教师：*学 院：理学院2011年5月摘 要matlab是 mathworks公司推出的一套科学计算软件，matlab的意思是矩阵实验室。matlab具有起点低、功能强大、易学易用以及兼有数值运算和符号运算功能的优点。利用matlab，绘图十分方便，它既可以绘制各种图形，包括二维图形和三维图形，还可以对图形进行修饰和控制。本文通过在matlab环境下编写通过科学计算解决经典物理问题，如力学、热学、电磁学中的一些常见问题。本文的思路主要是，先介绍经典物理习题，然后对习题进行分析，解答，再通过matlab软件进行编程，模拟实验结果。通过多次验证。得到所需答案。再通过图形绘制，形象的描绘出图形，与预期结果进行比较、验证。作出总结。本文展示的matlab软件在解决物理问题中的应用。关键词：力学；热学；电磁学；matlab程序abstract.mathworks matlab is introduced in a scientific computing software, matlab means matrix laboratory . matlab has a low starting point, powerful, easy to use, and both numerical calculation and symbolic operation advantages. using matlab, the drawing is very convenient, both to draw various graphics, including the two-dimensional graphics and three-dimensional graphics, graphics can also be modified and controlled. this article written by the matlab environment to solve by classical physics scientific computing problems, such as mechanical, thermal, electromagnetics some common problems. the main idea of this paper is to introduce classical physics problems, and then exercises to analyze, answer, and then programmed by matlab software to simulate the experimental results. through multiple authentication. get the answers you need. and through graphics rendering, the image depicts the graphics, compared with the expected results to verify. conclusion. this article presents the matlab software to solve the problem of physics.key words：mechanics;heat;electromagnetism,;matlab 目 录引言11 力学问题31.1质点运动学31.1.1已知质点的运动方程求其速度和加速度31.1.2已知质点的运动方程求质点的轨迹41.1.3考虑空气阻力的抛射体运动51.1.4已知加速度求速度、运动方程和轨迹71.2卢瑟福散射（rutherford scattering)研究92 热学问题122.1理想气体物态方程122.2范德瓦耳斯方程132.2.1范德瓦耳斯气体等温线132.2.2临界参数153电磁学问题163.1求电偶极子在其所在平面产生的电场中任一点p的电位163.2由电位的表示式计算电场并画出等电位线和电场方向173.3带电粒子在电磁场中的运动19结论21参考文献22致谢2322 引言近几十年来，计算机技术的广泛应用已经深入地影响到社会的各个方面，大大加快了社会的变革进程，计算机的应用离不开计算语言，而计算语言本身也处于不断的发展之中。matlab是matrix laboratory （矩阵实验室）的缩写，它自从1984年由美国mathworks 公司推出以来，经过不断改进和发展，现已经成为国际公认的优秀的工程应用开发环境。matlab是一种广泛应用于工程计算及数值分析领域的新型高级语言。它以矩阵作为数据操作的基本单位，使得矩阵运算变得非常简捷、方便、高效。matlab提供了十分丰富的数值计算函数，而且matlab和著名的符号计算语言maple相结合，使得matlab具有符号计算功能。matlab的绘图功能也很强，它既可以绘制各种二维、三维图形，还可以对图形进行修饰和控制，以增强图形的表现效果。matlab具有编程语言的基本特征，使用matlab也可以像使用basic、fortran、c等传统编程语言一样，进行程序设计，而且简单易学、编程效率高。matlab包含基本部分和各种可选的工具箱，其基本部分构成了matlab的核心内容，而matlab工具箱扩充了其功能。应用范围也越来越广。物理模型的建立及其数学处理在物理学的教学中占有重要地位，而matlab在这方面具有独特的优势。因此，利用matlab这一先进的科学计算语言来辅助物理学的教学工作必将大大提高教学效率。另外，matlab起点低、功能强、易学易用以及兼有数值运算和符号运算功能的优点，可以初步掌握这门科学计算语言，并在整个物理学习过程中不断反复使用是完全必要和可行的。运动学的任务是描述随时间的推移物体空间位置的变动，不涉及物体间相互作用与运动的关系。本文在力学中主要讨论如何使用matlab描述质点理想模型的运动，最后引入伽利略变换，它和物理学一条基本原理即相对性原理密切相关。质点平面运动指质点在平面上的曲线运动。这时，质点经常改变运动方向，速度、加速度等物理量的矢量性更突出。如何选择坐标系的问题更加重要。本文在质点运动方面，主要讨论抛体运动，在理想情况下，受空气阻力、斜抛等得运动轨迹如何在matlab中体现出来。以及，已知速度、如何求加速度等。本文在热学方面主要处理了理想气体物态方程、范德瓦耳斯方程如何用matlab描述出来。理想气体，只要在足够宽广的温度、压强变化范围内进行比较精细的研究，就可发现，气体的物态方程相当复杂，而且不同气体所遵循的规律也有所不同。但在压强趋于零，其温度不太高也不太低的情况下，不同种类气体在物态方程上的差异可趋于消失，气体所遵从的规律也趋于简单。这种压强趋于零的极限状态下的气体称为理想气体。荷兰物理学家范德瓦耳斯在克劳修斯的论文的启发下，对理想气体的两条基本假定即忽略分子固有体积、忽略除碰撞外分子间相互作用力作出了两条重要修正，得出了能描述真实气体行为的范德瓦耳斯方程。在发现电现象2000多年之后，人们才开始对电现象进行定量的研究。1785年，库伦通过扭秤实验总结出两个静止电荷之间电相互作用的定量规律，通常称之为库仑定律。实验表明，静电力具有叠加性。原则上，库仑定律加上静电力的叠加原理可以求解任意带电体之间的静电力。实验也指出，试探电荷在场中所受的静电力与试探电荷电量之比反映了电场本身的性质，该比值被称为电场强度。电场强度也具有叠加性，由场强的定义加上场的叠加原理可以求解任意带电体的场强分布。本文在电磁学中，主要研究如何用matlab求解电偶极子，带电粒子在电场中运动的问题。本文在物理题目的选取上，主要是普遍、常见的问题，意在将计算语言和物理课程的学习结合起来，起到相辅相成的作用。在程序的编写中，也力求简洁。1 力学问题1.1质点运动学在一些问题中，若物体的形状和大小可以忽略，则可以把该物体视为具有一定质量的几何点，这就是所谓的质点。质点运动学的基本问题是；已知质点的运动学方程求质点的轨迹、速度和加速度；已知质点的速度或加速度求质点的运动方程和轨迹。下面，结合大家熟悉的几个具体例子来说明如何用matlab处理上述问题1.1.1已知质点的运动方程求其速度和加速度例：某质点的运动学方程为（单位：m，s），求t=1s时质点的速度矢量。解题分析质点的位置矢量为，质点的速度矢量为，质点速度大小和方向余弦分别为程序syms tr=-10,15*t,5*t3; %用数组表示位置矢量v=diff(r,t); %求速度v=sqrt(sum(v.2) %求速度矢量长度，即速度矢量的大小alpha=acos(v(1)/v); beta=acos(v(2)/v); gamma=acos(v(3)/v);%求速度矢量的方向角v1=subs(v,t,1), alpha=subs(alpha,t,1),beta=subs(beta,t,1),gamma=subs(gamma,t,1)%求t=1s时质点的速率和方向角，使用了置换命令的函数subs运行结果：v1=21.2132alpha=1.5708beta=0.7854gamma=0.78541.1.2已知质点的运动方程求质点的轨迹例：设一物体以抛射角，速度抛出，落点与射点在同一水平面，且不计空气阻力。求物体在空气中飞行的时间、落点距离和飞行的最大高度。解题分析：质点运动学，有 解出t，它就是落点时间.有两个解，只取其中的一个有效解，然后求最大飞行距离。matlab程序：clear all y0=0;x0=5; %取初始位置，为了画出竖抛运动，未将x0取在原点。v0=input(v0=);theta=input(theta=); %输入抛射速率和岀射角度v0x=v0*cosd(theta);v0y=v0*sind(theta); %输入初速度的x分量和y分量ay=-9.81;ax=0; %加速度的y分量和x分量tf=roots(ay/2,v0y,y0); %解出方程的根，求飞行时间。有两个解，只取有效解tf=max(tf); %落点时间t=0:0.1:tf; %为了画图，取时间数组y=y0+v0y*t+ay*t.2/2;x=x0+v0x*t+ax*t.2/2; %t时刻，质点的位置xf=max(x), %飞行达到的最远距离，即射程yf=max(y), %飞行中达到的最大高度grid on, hold on %画网格，保持图形plot(x,y), %画图，xlabel(x),ylabel(y) %坐标标注hold off仿真结果与讨论：运行该程序，例如，取初速度=30，岀射角分别取35,45,55,65，75,85，90，则可画出图1.1所示曲线，并在命令窗口中给出相应的射程、飞行时间和最大高度。图 1.1抛体的运动轨迹在上述程序中，我们设定了目标高度为零。我们还可对上述程序进行修改，使其能预先设定目标高度。1.1.3考虑空气阻力的抛射体运动例：一物体在有阻力的空气中作抛体运动，设抛体质量为m，初速度为（可设），受到的空气阻力大小与速率v的一次方成正比，b是比例系数。求抛体的运动轨迹和速度的x、y分量以及速率v随时间的变化。解题分析：以地面为参考系，以抛出点为原点建立直角坐标系。质点受重力和空气阻力的作用，其运动微分方程为令y(1)=x, y(2)=dx/dt, y(3)=y, y(4)=dy/dt, 将方程写成一阶微分方程组的形式再用命令函数ode45解此常微分方程组。matlab程序m=1; %为简单起见，取m=1.b=input(b=); %输入b值，例如，b=0.3.t,y=ode45(ex1,0:0.01:5,0,5,0,19,b,m);%使用了数值法解微分方程的命令函数ode45v=sqrt(y(:,2).2+y(:,4).2);subplot(2,1,1) %绘制子图plot(y(:,1),y(:,3) %绘制运动轨迹，即x-y曲线，注意：y(1)=x,y(3)=y.subplot(2,1,2)plot(t,y(:,2),t,y(:,4),t,v) %绘制速度的x，y分量以及速率v时间t的变化曲线。函数文件function ydot=ex1(t,y,flag,b,m)ydot=y(2); -b./m.*y(2).*sqrt(y(2).2+y(4).2); y(4);-9.8-b./m.*y(4).*sqrt(y(2).2+y(4).2);运行结果如图1.2所示。可以改变b值（例如b分别取0.3和0）来观察运动轨迹和速度的x分量、y分量 及速率v随时间的变化。（a)b=0.3(b)b=0图1.2有空气阻力时抛射体的运动轨迹及速度随时间的变化1.1.4已知加速度求速度、运动方程和轨迹例：一质点平面运动的加速度为。初始条件为。求质点轨迹。解题分析：将加速度对时间求积分可得速度，将速度对时间求积分可得位置坐标。在得到参数方程后，给定时间数组，就可以画出运动轨迹。matlab程序clear allsyms t a b v0x v0y x0 y0 vx vy ax ay t t1 t2;v0x=0;v0y=b; x0=a; y0=0; %初始条件ax=-a*cos(t); ay=-b*sin(t); %加速度的x分量和y分量。vx=v0x+int(ax,t,0,t1),vy=v0y+int(ay,t,0,t1), %速度的x分量和y分量x=a+int(vx,t1,0,t2); y=int(vy,t1,0,t2); %求参数方程x=vpa(subs(x,a,b,1,0.5),3)%使用了vpa计算数值；使用subs用数据替换a和b。y=vpa(subs(y,a,b,1,0.5),3)运行结果：vx=-a*sin(t1)vy=-b*cos(t1)x=cos(t2)y=500*sin(t2)%下面 绘制质点的轨迹 cleart2=0:0.1:2*pi;x=cos(t2);y=.500*sin(t2);plot(x,y),xlabel(x),ylabel(y)从图1.3给出了程序的运行结果，可知质点的运动轨迹是椭圆。图1.3：质点的轨迹1.2卢瑟福散射（rutherford scattering)研究例：卢瑟福等人发现用粒子轰击金箔时有些入射粒子散射偏转角很大，甚至超过。卢瑟福于1911年提出原子必有以带正电的核心，即原子核；此即原子结构的行星模型。已知粒子的质量为，以速度接近电荷为ze的重原子核，瞄准距离为b，如图所示。设原子核质量比粒子大很多，可以近似看作静止。1 求粒子接近重原子核最近距离d。2 画出粒子在不同初始条件下的轨道，并通过改变初始条件来研究影响散射角的因素。解题分析1 粒子受静电力始终指向重核中心，粒子在一平面内运动。设z轴垂直于此平面且通过重核中心，则粒子所受静电力对z轴的力矩为零，即对z轴的角动量守恒。粒子以速度运动，对z轴的角动量是，粒子最接近重核（距离为d）时，速度应与粒子至核的连线垂直，角动量为。于是或（1）在散射过程中，只有库仑斥力作用，故能量守恒。（2）其中，左边第二项是库仑斥力势能。联解（1）、（2）式，可得d的表达式。2 选择在直角坐标系，原点位于力心重核处。根据牛顿运动定律，粒子的运动方程在直角坐标中的投影方程为令，则上述方程组可写为令粒子沿ox方向入射，入射速率为，初始条件为。为了能得到多粒子的运动轨迹，程序中采用input函数给出不同初始条件。matlab程序1 求粒子接近重原子核最近距离d。syms v v0 b k z e m d;d,v=solve(v=v0*b/d,m*v2/2+k*z*e2/d=m*v02/2,d,v)运行结果：d=1/2/m/v02*(2*k*z*e2+2*(k2*z2*e4+m2*v04*b2)(1/2)v=2*v03*b*m/(2*k*z*e2+2*(k2*z2*e4+m2*v04*b2)(1/2)即2 画出粒子在不同初始条件下的轨道y0=input(请输入初始条件:); %例如，可输入：-20 1 10 0;line(0,0,marker,.,markersize,50,color,r);text(2,0,靶粒子);hold ont,y=ode23(ex2f,0:.1:42,y0,3);axis(-20 20,-20 20)plot(y(:,1),y(:,3),hold on以下是独立的函数文件，文件名为ex2f.m,其中。function ydot=ex2f(t,y,flag,p)ydot=y(2);p*y(1)./sqrt(y(1).*y(1)+y(3).*y(3).3;y(4);p*y(3)./sqrt(y(1).*y(1)+y(3).*y(3).3;运行结果如图1.4所示。图1.4：粒子的散射轨道2 热学问题2.1理想气体物态方程理想气体是将实际气体外推到压强趋于零的极限情况下得到的一个理想模型。1857年，克劳修斯进一步提出了理想气体的微观模型，并通过计算气体的压强得到了理想气体的物态方程。而在此之前，理想气体物态方程是由气体三大实验定律外推得到的。例：编写一个绘制带有等高线的理想气体状态方程pv=rt的曲面。解题分析理想气体的物态方程为。其中，p，v，分别为气体的压强、体积和摩尔数，r为气体普适常数，其值为r=8.31j/mol*k.图2.1：理想气体状态方程曲面图matlab程序clearr=8.31;p=(1:20).*1e5;v=(1:20)*1e-3;v,p=meshgrid(v,p);t=p.*v./r;meshc(v,p,t),xlabel(v),ylabel(p),zlabel(t),运行结果如图2.1所示。实验指出，理想气体状态方程在一定程度上反映了真实气体的性质，但对低温和高密度状态下的气体以及气体和液体之间的相变却无能为力，因而是一个理想的“永久气体”状态方程。2.2范德瓦耳斯方程范德瓦耳斯方程1873年，荷兰物理学家范德瓦耳斯（van der waals)在克劳修斯的理想气体模型和安德鲁斯发现的临界点现象的启发下，考虑了分子体积和分子间吸引力这两个因素，对理想气体进行了修正，得到了能描述真实气体行为的范德瓦耳斯方程：其中，常数a和b分别是1mol范氏气体的压强修正系数和体积修正系数，其数值随气体种类的不同而异。下表1列出了几种气体的a，b值及临界参量。 表1：几种气体的a、b值及临界参量气体 atc(k)ar0.13490.0332747.9789150.86760.09680.36060.0428072.9079304.06400.1284co0.14560.0395434.4925132.89470.1186he0.0034140.023712.24925.19650.07110.024610.0266812.804933.28960.0800ne0.020910.0169726.892244.4680.05090.13500.0386433.4885126.08970.11590.13630.0318449.7950154.49200.09552.2.1范德瓦耳斯气体等温线例:编写一个绘制范德瓦尔斯气体等温线的程序，要求输入温度值后便可画出相应的等温线。解题分析以二氧化碳为例，从表1查得，由范德瓦尔斯方程可绘制等温线簇，温度选取如图所示。matlab程序v=(0.06:0.001:1).*1e-3;t=input(t=);b=0.0428e-3;a=0.3606;r=8.31;p=r.*t./(v-b)-a./v.2;grid on,plot(v,p),axis(0,0.4e-3,-2e7,6e7)hold on;运行结果如图2.2所示图2.2：范德瓦耳斯气体等温线范德瓦尔斯方程不仅对气体性质的描述优于理想气体物态方程，而且还能描述液相及气、液两相转变的性质以及临界点的特征。2.2.2临界参数范德瓦尔斯等温线中有一个特殊的状态临界点。临界点所对应的压强、体积和温度分别称为临界压强、临界体积和临界温度。在临界点所发生的气液相变与在低于临界温度时的相变完全相同，属于二级相变；而低于临界点是的气液相变属于一级相变。在临界点以上，气体是不能够通过等温压缩被转变为液相的。系统在临界点具有许多特殊性质，称为临界现象。下面来介绍临界点的确定。例:由范德瓦尔斯物态方程求临界参量、。解题分析从图可以看出，临界点是一拐点，它同时满足下列条件：，利用上述拐点条件，将范德瓦尔斯方程对求导并联解方程，便可求得范德瓦尔斯气体的三个临界参量。matlab程序clearsyms a b r td1=diff(p+a/v2)*(v-b)-r*t,v);d2=diff(d1,v);pc,vc=solve(d1,d2,v,p)tc=solve(pc+a/vc2)*(vc-b)-r*t,t)运行结果pc=1/27*a/b2vc=3*btc=8/27*a/b/r即；，3电磁学问题3.1求电偶极子在其所在平面产生的电场中任一点p的电位例：已知电偶极子中两电荷-q和+q的距离为。计算中可取。解题分析设场点p到的距离为和，则单独存在时p点的电位分别为由电位叠加原理，电偶极子产生的电场在p点的电位为matlab程序clear;q=1.6e-19; %单位电荷电量c0=9e9;l=3.0; %偶极子正负电荷之间的距离lx=-5:0.5:5;y=x;x,y=meshgrid(x,y);r1=sqrt(x-1./2).2+(y-0).2); %电荷距离空间p（x，y）点的距离r2=sqrt(x+1./2).2+(y-0).2);u=q.*c0.*(1./r1-1./r2); %求出空间任意一点p（x，y）的电位plot(-1/2,0,ro,-1/2,0,r-) %标出负电荷hold on,plot(1/2,0,ro,1/2,0,r+) %标出正电荷c=contour(x,y,u,k-);clabel(c); %画等位线并标出电位值axis(square)运行结果如图3.1所示 3.1电位梯度与电场强度电位是标量，它在空间中每点都有一定的数值，所以电位是标量场。标量场在空间中沿不同方向的变化率称为梯度，对电位场而言称为电位梯度，用grad u或来表示。可以证明，电位梯度和电场强度e的关系为利用上式，可从已知的电位分布求电场强度。3.2由电位的表示式计算电场并画出等电位线和电场方向解题分析如果已知空间的电位分布则空间的电场强度为按照本题的要求，可利用读入字符串的指令input(u(x,y)=,s)来输入电位方程。在matlab中，梯度函数的调用格式为gradient（），它是靠数值微分得到的。因此，空间观测点应取得密一些，以获得较高的精度。matlab程序clear allu=input(请输入电位方程，u=(x,y)=,s); %例如，取u(x,y)=log(x.2+y.2)。xmax=5;ymax=5;ngrid=20; %绘图区从x=-xmax到xmax，网格线数为20xplot=linspace(-xmax,xmax,ngrid); %绘图用x的数组x,y=meshgrid(xplot); %x，y取同样范围，生成二维网格uplot=eval(u); %执行输入的字符串u,计算各点u的值explot,eyplot=gradient(-uplot); %电场等于电位的负梯度clf;subplot(1,2,1),meshc(uplot); %划分子图；绘制含等位线的三维曲面xlabel(x);ylabel(y);zlabel(u);subplot(1,2,2),axis(-xmax,xmax,-ymax,ymax); %规定等位线的范围cs=contour(x,y,uplot); %画等位线，cs是等位线值clabel(cs); %标出等位线的值hold on;quiver(x,y,explot,eyplot); %保持图形，在原图形上叠加矢量场图xlabel(x);ylabel(y);hold off;图3.2：的电位分布与电场分布运行上述程序，所得结果如图3.2所示。3.3带电粒子在电磁场中的运动例：设质量为m，带电量为q的粒子在电磁场中的运动微分方程为选场中某点为原点，以为方向，沿方向，建立坐标系。令，上式的投影方程为令，上述方程可改写为下列一阶微分方程组：matlab程序%符号法求离子运动微分方程的特解并绘图clearsyms w x y z t b e m q;e=input(e=);b=input(b=); %输入e和b值x,y,z=dsolve(d2x=q*b/m*dy,d2y=q*e/m-q*b/m*dx,d2z=0,x(0)=0,y(0)=0,z(0)=0,dx(0)=0.01,dy(0)=6,dz(0)=0.01);%初始条件取x(0)=y(0)=z(0)=0,dx(0)=0.01,dy(0)=0.01q=1.6e-2;m=0.02x=subs(x,y,z);x=x(1),y=x(2),z=x(3),ezplot3(x(1),x(2),x(3)运行上述程序，例如，取e=4，b=8可得下列特解并给出图3.3x =-15/16*cos(32/5*t)-49/640*sin(32/5*t)+1/2*t+15/16y =15/16*sin(32/5*t)-49/640*cos(32/5*t)+49/640z =1/100*t(a)e=4,b=8(b)e=0.01,b=8(c)e=8,b=1图3.3：带电粒子在电磁场中的运动下面我们给出一段用数值方法求解该问题的程序，用以比较。q=1.6e-2;m=0.02;b=2;2;0;e=1;0;1;figurestrd1=eneq 0,bneq 0;strd2=e=0,bneq 0;strd3=eneq 0,b=0;for i=1:3t,y=ode23(ex3f,0:0.1:20,0,0.01,0,6,0,0.01,q,m,b(i),e(i);axes(unit,normalized