Matlab 在电磁场中的应用ppt课件.ppt

Matlab在电磁场中的应用 2009 10 8 Matlab是美国Mathworks公司于80年代推出的大型数学软件 通过多年的升级换代 现在已发展成为集数值计算 符号计算 可视化功能以及诸多的工具箱为一体的大型科学计算软件 它已广泛应用于科研院所 工程技术等各个部门 并成为大学生 研究生必备的工具软件 Matlab语言以矩阵为基本运算单元 界面极其友好 初学者只要具备必要的数学基础就可以以令人难以置信的速度将其学到手 它对许多物理学的典型问题的计算和处理只需很少的几条语句就能做到 特别是它强大的作图及动画功能能够准确而科学地反应物理学问题的本质 语句也是简洁清晰 Matlab特点 功能强大的数值运算功能 强大的图形处理能力高级但简单的程序环境丰富的工具箱与模块集易于扩充 预备知识 一MATLAB运算1 算术运算 1 基本算术运算MATLAB的基本算术运算有 加 减 乘 右除 左除 乘方 注意 运算是在矩阵意义下进行的 单个数据的算术运算只是一种特例 2 点运算在MATLAB中 有一种特殊的运算 因为其运算符是在有关算术运算符前面加点 所以叫点运算 点运算符有 和 两矩阵进行点运算是指它们的对应元素进行相关运算 要求两矩阵的维参数相同 2 关系运算MATLAB提供了6种关系运算符 大于 大于或等于 等于 不等于 它们的含义不难理解 但要注意其书写方法与数学中的不等式符号不尽相同 3 逻辑运算MATLAB提供了3种逻辑运算符 与 或 和 非 二程序控制结构 2选择结构 if语句在MATLAB中 if语句有3种格式 1 单分支if语句 if条件语句组end当条件成立时 则执行语句组 执行完之后继续执行if语句的后继语句 若条件不成立 则直接执行if语句的后继语句 双分支if语句 if条件语句组1else语句组2end当条件成立时 执行语句组1 否则执行语句组2 语句组1或语句组2执行后 再执行if语句的后继语句 例计算分段函数的值 程序如下 x input 请输入x的值 ifx 0y x sqrt pi exp 2 elsey log x sqrt 1 x x 2 endY 多分支if语句 if条件1语句组1elseif条件2语句组2 elseif条件m语句组melse语句组nend语句用于实现多分支选择结构 3循环结构for语句for语句的格式为 for循环变量 表达式1 表达式2 表达式3循环体语句end其中表达式1的值为循环变量的初值 表达式2的值为步长 表达式3的值为循环变量的终值 步长为1时 表达式2可以省略 while语句while语句的一般格式为 while 条件 循环体语句end其执行过程为 若条件成立 则执行循环体语句 执行后再判断条件是否成立 如果不成立则跳出循环 三几个绘图命令1 plot函数 用来绘制线形图形plot y 当y是实向量时 以该向量元素的下标为横坐标 元素值为纵坐标画出一条连续曲线 这实际上是绘制折线图 plot x y 其中x和y为长度相同的向量 分别用于存储x坐标和y坐标数据 plot x y s contour函数 用来绘制等高线图形ezplot函数 对于显式函数f f x 在默认范围 2 x 2 上绘制函数f x 的图形 对于隐式函数f f x y 在默认的平面区域 2 x 2 2 y 2 上绘制函数f x y 的图形 2 具有两个纵坐标标度的图形在MATLAB中 如果需要绘制出具有不同纵坐标标度的两个图形 可以使用plotyy绘图函数 调用格式为 plotyy x1 y1 x2 y2 其中x1 y1对应一条曲线 x2 y2对应另一条曲线 横坐标的标度相同 纵坐标有两个 左纵坐标用于x1 y1数据对 右纵坐标用于x2 y2数据对 3三维曲线plot3函数与plot函数用法十分相似 其调用格式为 plot3 x1 y1 z1 选项1 x2 y2 z2 选项2 xn yn zn 选项n 其中每一组x y z组成一组曲线的坐标参数 选项的定义和plot函数相同 当x y z是同维向量时 则x y z对应元素构成一条三维曲线 当x y z是同维矩阵时 则以x y z对应列元素绘制三维曲线 曲线条数等于矩阵列数 一 单电荷的场分布 单电荷的外部电位计算公式等位线就是连接距离电荷等距离的点 在图上表示就是一圈一圈的圆 而电力线就是由点向外辐射的线 比较简单 这里就不再赘述 theta 0 0 01 2 pi r 0 10 x sin theta r y cos theta r plot x y b x linspace 5 5 100 fortheta pi 40pi 4 y x tan theta holdon plot x y endgridon 单电荷的等位线和电力线分布图 二 点电荷电场线的图像 考虑一个三点电荷系所构成的系统 如图所示 其中一个点电荷 q位于坐标原点 另一个 q位于y轴上的点 最后一个 2q位于y轴的 点 则在xoy平面内 电场强度应满足 y x q q 2q 任意条电场线应满足方程 1 求解式 1 可得 2 这就是电场线满足的方程 常数C取不同值将得到不同的电场线 解出y f x 的表达式再作图是不可能的 用Matlab语言即能轻松的做到这一点 如图2所示 其语句是 symsxy 设置x y变量 forC 0 0 1 3 0ezplot 2 y 1 sqrt y 1 2 x 2 y sqrt y 2 x 2 y 1 sqrt y 1 2 x 2 C 5 5 0 1 holdon end其中取了a 1 C 0 0 1 0 2 3 0 三 线电荷产生的电位 设电荷均匀分布在从z L到z L 通过原点的线段上 其密度为q 单位C m 求在xy平面上的电位分布 点电荷产生的电位可表示为是一个标量 其中r为电荷到测量点的距离 线电荷所产生的电位可用积分或叠加的方法来求 为此把线电荷分为N段 每段长为dL 每段上电荷为q dL 看作集中在中点的点电荷 它产生的电位为然后对全部电荷求和即可 把xy平面分成网格 因为xy平面上的电位仅取决于离原点的垂直距离R 所以可以省略一维 只取R为自变量 把R从0到10米分成Nr 1点 对每一点计算其电位 matlab程序clearall L input 线电荷长度L N input 分段数N Nr input 分段数Nr q input 电荷密度q E0 8 85e 12 C0 1 4 pi E0 L0 linspace L L N 1 L1 L0 1 N L2 L0 2 N 1 Lm L1 L2 2 dL 2 L N R linspace 0 10 Nr 1 fork 1 Nr 1Rk sqrt Lm 2 R k 2 Vk C0 dL q Rk V k sum Vk end max V min V plot R V grad 输入 线电荷长度L 5分段数N 50分段数Nr 50电荷密度q 1可得最大值和最小值为 ans 1 0e 010 9 31990 8654 线电荷产生的静电位分布图 四 计算平面上N个电荷之间的库仑引力 建模 由库仑定律 其分量的公式可以写成 编写程序时 先输入电荷的数目 各电荷的坐标及电荷量 再选一个电荷 求其它电荷对它的作用力 叠加求合力 再选下一个电荷 依次类推 Matlab程序 clearall N input 输入电荷数目N foric 1 N 输入给定条件fprintf n对电荷 g n ic rc input 输入电荷位置 x y 米 x ic rc 1 电荷ic的x坐标y ic rc 2 电荷ic的y坐标q ic input 输入电荷量 库仑 end E0 8 85e 12 真空中的常数C0 1 4 pi E0 合并常数foric 1 N 循环计每个电荷所受的力Fx 0 0 Fy 0 0 forjc 1 Nif ic jc xij x ic x jc yij y ic y jc Rij sqrt xij 2 yij 2 Fx Fx C0 q ic q jc xij Rij 3 Fy Fy C0 q ic q jc yij Rij 3 endend fprintf 其它电荷作用在电荷 g上的合力为 n ic fprintf x 分量 gN n Fx fprintf y 分量 gN n Fy end本程序注意学会循环提示并输入参数的方法 以及用双循环解决较复杂的计算过程的编程问题 输入已知条件 输入电荷数目N 3 对电荷 1输入电荷位置 x y m 12 输入电荷量 库仑 2 对电荷 2输入电荷位置 x y m 11 输入电荷量 库仑 1 对电荷 3输入电荷位置 x y m 33 输入电荷量 库仑 3 计算结果 其它电荷作用在 1上的合力为 X 分量为 9 65102e 009NY 分量为1 31581e 010其它电荷作用在 2上的合力为 X 分量为 2 38431e 009NY 分量为 2 03679e 010其它电荷作用在 3上的合力为 X 分量为 1 20353e 010NY 分量为7 20982e 009 有限差分法处理电磁场问题 在很多实际情况下 往往我们不知道电荷的分布情况 而只知道边界的电位 例如在静电场中的导体 边界是等位面 但是表面上的电荷分布往往不一样 而且很难求得 如果我们想求导体附近的电场 这里就介绍用差分法解电场的方法 解决这个问题的关键是对电位使用以下结论 在一个没有电荷的区域 给出一个点的电位等于周围点的电位数值的平均值 我们使用高斯定律以及以下公式来证明这个结论 根据电压V的电场分量 T8 1 我们将集中讨论该情况 其中电位只取决于两个坐标 x和y 一个例子是一个带电的长圆柱体 在一个点的电势只依赖于这一点在平面垂直于圆柱体的轴线 而不是z坐标 对于这样一个二维的情况下 考虑一个点P的坐标 x Y Z 并在一个由高斯表面封闭的立方体的一面长度是2 l 中心在P 图T8 1 如果立方体的内部没有电荷 通过立方体的电通量 E等于零 由方程 T8 1 可知Z轴的电场分量为零 因为电势V并不是z的函数 因此 并没有通过高斯表面的平行于xy 平面的电通量 由于是一个小的立方体 通过其他四个面每通量有一个良好的逼近 等于在每面的中心和每一面的 2 l 2的E的垂直分量的乘积 总流量 等于0 可以表示为 T8 2 T8 1 在一个没有电荷的区域 点P的电位数值等于P点周围电位值的平均值 使用 T8 1 我们可以近似的写出电场的各分量 我们可以得到P点的电位是 T8 3 T8 4 总之 P点的电位值等于P点周围点的电位值的平均值 前提是非常小 下面是用计算机绘图的过程 1 选取一个正电压 2 选取横纵线上网格点数m 则在内场区的网格点数为3 V j k 为网格上点的坐标 j从1到m k从1到m 4 开始循环k从1到m 5 开始循环j从1到m 6 给网格上的点赋初值 k 1时 V j k j 1 j m或者k m时V j k 0 其余的V j k 0 7 j和k都结束循环 8 设置一个值 为最小精度 这里设为cha 0 01 9 设置一个变量 初值为0 这里设为delta 0 10 k开始循环 从2到们m 1 11 j开始循环 从2到们m 1 就是内场的点 不包含边界 12 设置一个变量使之等于 13 设置一个变量d 使之等于如果d delta 则把d的值赋给delta 14 把的值赋给以前的V j k 15 结束j k的循环 16 如果delta的值大于前面设置的精度cha 则至少有一个网格上的点两次计算之间的差值大于所设的精度 则程序从第10步开始重新计算 直到delta的值小于cha为止 17 输出网格上点的电压值 18 程序结束 计算场点电压值程序如下 m 12fork 1 mforj 1 mifk 1V j k 1 elseif j 1 j m k m V j k 0 elseV j k 0 5 endendend cha 0 01 delta 0 n 0 while 1 n n 1 fork 2 m 1forj 2 m 1Vnew j k 1 4 V j 1 k V j 1 k V j k 1 V j k 1 d abs Vnew j k V j k V j k ifd deltadelta d endV j k Vnew j k endendifdelta100 break enddelta 0 end 程序计算出的结果如下 代入m 22 绘图程序 k 1 m j 1 m DX DY gradient V 4 4 holdonquiver k j DX DY 2 holdoffk 1 m j 1 m DX DY gradient V 4 4 A DX 2 DY 2 0 5 DA DB gradient A 4 4 holdonquiver k j DA DB 2 holdoff 电场图 电力线 通过以上的初步探讨 可以看出 在有限差分法的计算中 或者说是在工程计算和数据处理中 确实具备如下优点 1 较其它高级程序设计语言 Matlab程序语言的规则更为接近数学表示 2 语句简洁明了 表意却出乎意料的丰富 比如 上面程序中函数mesh contour