计算物理2013

1、计算物理基础计算物理基础 九江学院理学院 2013.22022-5-9计算物理基础2计算物理是以电子计算机为工具、采用数学方法解决物理问题的应用科学。本课程的目的在于对计算物理进行一些入门指导，使大家在学完本课程后，能利用Matlab软件进行物理问题求解，以及分析处理其它可以数值求解的问题。绪论绪论课程目的2022-5-9计算物理基础3掌握计算物理的概念和方法；掌握几类计算方法的基础或基本原理；了解这些方法在若干物理学分支中的具体应用。计算物理的实践性非常强，上机是本课程的一个有机组成部分本课程需具备高等数学和线性代数基本知识课程要求2022-5-9计算物理基础4第一章第一章 MatlabMa

2、tlab简介简介第二章第二章 迭代与分形图形迭代与分形图形第三章第三章 数值微分与数值积分数值微分与数值积分第四章第四章 数据处理数据处理第五章第五章 解常微分方程解常微分方程第七章第七章 解偏微分方程解偏微分方程第八章第八章 蒙特卡罗方法蒙特卡罗方法 目目 录录2022-5-9计算物理基础50.10.1、什么是计算物理？、什么是计算物理？0.20.2、计算物理、计算物理 的起源、形成与发展的起源、形成与发展0.30.3、计算物理的进一步发展、计算物理的进一步发展 从计算物理到科学计算、战略计算从计算物理到科学计算、战略计算0.40.4、计算物理的工作流程、计算物理的工作流程第一章第一章 绪论

3、 MatlabMatlab简介简介2022-5-9计算物理基础61 什么是计算物理？物理学有几大门类？传统物理学分为理论物理与试验物理两大分支理论物理实验物理计算物理？2022-5-9计算物理基础7 理论物理是分析的科学，它从一系列的基本理论物理是分析的科学，它从一系列的基本原理和基本假设出发，列出相应的数学方程，原理和基本假设出发，列出相应的数学方程，运用传统的或现在的数学方法求出问题的显式运用传统的或现在的数学方法求出问题的显式解析解，用这些解析解的结论去解释物理现象，解析解，用这些解析解的结论去解释物理现象，预见新的现象，指导实验。预见新的现象，指导实验。1 1 什么是计算物理？什么是计

4、算物理？2022-5-9计算物理基础8 实验物理是从实验观测出发，发现新的物理现象，为理论物理提供总结新的物理规律的素材，检验理论物理的假设或理论物理预言的正确程度和适用范围等1 什么是计算物理？2022-5-9计算物理基础9计算物理是伴随着电子计算机的出现和发展而逐步形成的一门新兴的边缘学科。是以电子计算机为工具、采用数学方法解决物理问题的应用科学。它是物理、数学和计算机三者相结合的产物。1 什么是计算物理？2022-5-9计算物理基础102 计算物理的起源、形成与发展 传统的物理学：理论物理，实验物理，都离不开数值计算，如海王星的发现及其轨道计算就是一个典型例子。 但早期的计算仅使用人力或

5、简单的计算工具，其功能和效率都极其有限。这种计算不能成为一个学科分支。 牛顿力学方程只有二体问题是可解得，三体以上的问题折磨了全世界许多优秀的数学家和理论物理学家，仍然没有解析解。 量子力学的薛定谔方程，除了氢原子和简谐振子外没有一个真实的物理问题可以找到解析解。2022-5-9计算物理基础11 20世纪40年代初，在由于战争的需要开始了核武器研制。涉及的问题：流体动力学过程、核反应过程、中子输运过程、光辐射输运过程、物态变化过程等；都是十分复杂的非线性方程组，不可能用传统的解析方法求解。 由于需要在短时间内进行大量复杂的数值计算，从而促使了计算机的延生和新物理学科的形成。2 计算物理的起源、

6、形成与发展2022-5-9计算物理基础12 科学家们从原子弹设计中使用计算机求解复杂物理问题取得成功而得到启示，迅速将这种方法推广应用到物理学的其他领域：天体物理、大气物理、等离子体物理、核物理、原子分子物理、固体物理、统计物理和基本粒子物理等，而且还应用到气象预报、水利、海洋、地震、石油、化工甚至人体科学等各个科学技术领域。2 计算物理的起源、形成与发展2022-5-9计算物理基础13 1963年，美国的Beini，Alder等人开始编辑出版计算物理方法丛书，内容涉及统计物理、量子力学、流体力学、核物理、天体物理、固体物理、等离子体物理、地球物理和大气环流等。 1966年，Journal o

7、f Computational Physics在美国创刊；1969年，Computer Physics Communication在西欧创刊。 1977年，美国和西欧的学者开始编辑出版计算物理施普林格系列丛书，到1988年已出17本；2 计算物理的起源、形成与发展2022-5-9计算物理基础141991年，以美国总统的名义提出“高性能计算与通信计划”。投资重点(43%)是发展先进的软件技术与并行算法，关键技术是可扩展的大规模并行计算。1993年美国总统发布“发展信息高速公路”的总统令1994年美国总统发布“建立国家（地球）空间数据基础设施”的总统令。所有这些计划，都是为大规模科学计算创造条件，

8、促使科学计算高速发展。3 3 计算物理的进一步发展计算物理的进一步发展 从计算物理到科学计算、战略计算从计算物理到科学计算、战略计算2022-5-9计算物理基础154 计算物理的特征计算物理的研究内容（计算机实验） 凡是局部瞬时的物理规律已知或被假设，要想求得大范围长时间的物理现象的发展过程，便属于计算物理学的范围。从局部关系到大范围依赖于计算机的大容量由瞬时规律发展为长时间的过程依赖于计算机的高速度。2022-5-9计算物理基础16计算物理相对于理论物理的优越性理论物理中利用数学方程组求解物理问题时，通常将问题大加简化，这些简化包括：复杂问题只考虑少数主要因素：质点，黑体近似等动态过程只考虑

9、最后达到的静态状况：热平衡等将非线性因素硬作线性化处理将变系数硬作常系数处理将复杂的边界简化为规则的边界等等4 计算物理的特征2022-5-9计算物理基础17计算物理提出计算预测给出模拟结果提供计算数据提供模拟结果检验计算预测提供方程实验物理提供实验数据解释结果理论物理检验理论预测提供实验数据提出理论预测给出理论解释计算物理与传统物理的联系4 计算物理的特征实验结果实验结果计算机计算结果计算机计算结果2022-5-9计算物理基础19计算物理方法区别于计算数学方法的特点：1）计算物理从物理问题出发，以物理结论为结果，以与实验数据的对比为其结束；而计算数学则是从数学方程出发，以求得方程的近似解告终

10、。计算物理工作者选用计算方法时要考虑算法和结果的物理意义；而计算数学工作者最感兴趣的是算法的逼近阶，计算精度和稳定性等问题。4 计算物理的特征国内超级计算中心发展掠影国家超级计算长沙中心（建设中，将完工）国家超级计算长沙中心（建设中，将完工）25Matlab 简介q Matlab: Matrix Laboratory 矩阵实验室q Matlab 的发展的发展l 1980年，Moler 教授用 Fortran 语言编写了集命令翻译、 科学计算于一身的一套交互式软件系统。l 1984年，Moler 等成立了 The MathWorks 的公司，用 C 语言完全改写 Matlab，并推出第一个商业版

11、。l 增添图形图像处理、符号运算、以及与其他流行软件 的接口功能，使得 Matlab 的功能越来越强大。l 到九十年代，在国际上 30 几个数学类科技应用软件中， Matlab 在数值计算方面独占鳌头。26q 目前，Matlab 已成为世界顶尖的数学应用软件，以其强大的工程计算、算法研究、工程绘图、应用程序开发、数据分析和动态仿真等功能，在航空航天、机械制造和工程建筑等领域发挥着越来越重要的作用。就影响而言，至今仍然没有一个别的计算软件可与 Matlab 匹敌。 q Matlab 的发行的发行1984年，Matlab 1.0 （DOS版，182K，20多个函数）1992年，Matlab 4.0

12、 （93年推出Windows版，加入 simulink）1994年，Matlab 4.2（得到广泛重视和应用）1999年，Matlab 5.3（真正实现32位运算）2002年，Matlab 6.5（采用JIT加速器）2004年，Matlab 7.0自2006年起，Matlab每年更新两次Matlab 简介27Matlab 的特点与功能q Matlab 具有很强的数值计算功能具有很强的数值计算功能l Matlab 以以矩阵矩阵作为数据操作的基本单位，作为数据操作的基本单位， 但无需预先指定矩阵维数（但无需预先指定矩阵维数（动态定维动态定维）l 按照按照 IEEE 的数值计算标准进行计算的数值计算

13、标准进行计算l 提供十分丰富的数值计算函数，方便计算，提高效率提供十分丰富的数值计算函数，方便计算，提高效率l Matlab 命令与数学中的符号、公式非常接近，命令与数学中的符号、公式非常接近， 可读性强，容易掌握可读性强，容易掌握q Matlab 是一个交互式软件系统是一个交互式软件系统输入一条命令，立即就可以得出该命令的结果输入一条命令，立即就可以得出该命令的结果28Matlab 的特点与功能q Matlab 符号计算功能符号计算功能Matlab 和著名的符号计算语言和著名的符号计算语言 Maple 相结合相结合q Matlab 的编程功能的编程功能Matlab具有具有程序结构控制程序结构

14、控制、函数调用函数调用、数据结构数据结构、输入输输入输出出、面向对象面向对象等程序语言特征，而且等程序语言特征，而且简单易学简单易学、编程效率编程效率高高。通过。通过 Matlab 进行编程完成特定的任务进行编程完成特定的任务q Matlab 的绘图功能的绘图功能Matlab提供丰富的绘图命令，提供丰富的绘图命令，很方便实现数据的可视化很方便实现数据的可视化29q Matlab 丰富的工具箱（丰富的工具箱（toolbox）Matlab 的特点与功能根据专门领域中的特殊需要而设计的各种可选工具箱根据专门领域中的特殊需要而设计的各种可选工具箱q Matlab 的的 Simulink 动态仿真集成环

15、境动态仿真集成环境提供建立系统模型、选择仿真参数和数值算法、启动仿提供建立系统模型、选择仿真参数和数值算法、启动仿真程序对该系统进行仿真、设置不同的输出方式来观察真程序对该系统进行仿真、设置不同的输出方式来观察仿真结果等功能仿真结果等功能Symbolic Math PDEOptimizationSignal processImage ProcessStatisticsControl SystemSystem Identification 31Matlab 的工作界面命令窗口命令窗口当前工当前工作目录作目录当前工当前工作空间作空间输入命令的输入命令的历史记录历史记录命令命令提示符提示符37各种

16、format 格式格式解释例format短格式（缺省显示格式），同short3.1416format short短格式（缺省显示格式），只显示5位3.1416format long长格式，双精度数15位，单精度数7位3.14159265358979format short e短格式e方式（科学计数格式）3.1416e+000format long e长格式e方式3.141592653589793e+000format short g短格式g方式3.1416format long g长格式g方式3.14159265358979format compact压缩格式format loose自由格式f

17、ormat + / format bank / format rat / format hex (详情查看联机帮助)38变量的存储q 存储当前工作空间中的变量存储当前工作空间中的变量u save 将所有变量存入文件将所有变量存入文件 matlab.matu save mydata 将所有变量存入将所有变量存入指定文件指定文件 mydata.matq 存储指定的变量存储指定的变量u save mydata.mat 将所有变量存入将所有变量存入文件文件 mydata.matsave 文件名文件名 变量名列表变量名列表例例： save mydata A x z 变量名列表中各变量之间用变量名列表中各

18、变量之间用空格空格分隔分隔39变量的读取q 将数据文件中的变量载入当前工作空间将数据文件中的变量载入当前工作空间u load mydata 载入数据文件中的所有变量载入数据文件中的所有变量u load mydata A x 从数据文件中提取指定变量从数据文件中提取指定变量q 清除当前工作空间中的变量清除当前工作空间中的变量u clear 清除当前工作空间中的所有变量清除当前工作空间中的所有变量u clear A x 清除指定的变量清除指定的变量48q 输出格式输出格式Matlab 的输出u Matlab 以双精度执行所有的运算，运算结果可以以双精度执行所有的运算，运算结果可以在在屏幕上输出屏幕

19、上输出，同时，同时赋给指定变量；赋给指定变量；若无指定变量，则系若无指定变量，则系统会自动将结果赋给变量统会自动将结果赋给变量 “ans” u Matlab 中数的输出格式可以通过中数的输出格式可以通过 format 命令指定命令指定format 只改变变量的输出格式，只改变变量的输出格式，但不会影响变量的值！但不会影响变量的值！57switch语句语句58（一）选择结构 选择结构的语句有if语句和switch语句。 1 if语句语句 格式一： if 条件 语句组 endn格式二： if 条件n 语句组1n elsen 语句组2n end59(一)选择结构格式三： if 条件1 语句组1 el

20、seif 条件2 语句组2 elseif 条件m 语句组m else 语句组m+1 end60【例4】 输入三角形的三条边，求面积。 A=input(请输入三角形的三条边：); if A(1)+A(2)A(3) & A(1)+A(3)A(2) & A(2)+A(3)A(1) p=(A(1)+A(2)+A(3)/2; s=sqrt(p*(p-A(1)*(p-A(2)*(p-A(3); disp(s); else disp(不能构成一个三角形。) end 运行： 请输入三角形的三条边：4 5 6 9.9216(一)选择结构61【例5】 输入一个字符，若为大写字母，则输出其后继字符，

21、若为小写字母，则输出其前导字符，若为其他字符则原样输出。 c=input(,s); if c=A & c=a& c No=input(Please input your choice! );Please input your choice! 1 switch Nocase 0disp(return to main menu);case 1disp(She is a girl);case 2disp(He is a boy);otherwisedisp(I cant determine) endShe is a girl65c) 循环结构循环结构while 语句语句:for语句语

22、句:for 变量=初值：增量：结束值 程序模块；end66 sum=0; i=0; while(i sumsum = 14196774举例举例： for I=1:10 A(I)=1/(I+1) ; end AA = Columns 1 through 7 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 Columns 8 through 10 0.1111 0.1000 0.090967break语句语句循环结构循环结构continue 语句语句循环结构循环结构d) 其它与流程控制有关的语句其它与流程控制有关的语句ifif68 3 循环的嵌套循

23、环的嵌套如果一个循环结构的循环体又包括一个循环结构，就称为循环的嵌套，或称为多重循环结构。多重循环的嵌套层数可以是任意的。可以按照嵌套层数，分别叫做二重循环、三重循环等。处于内部的循环叫作内循环，处于外部的循环叫作外循环。三程序设计(二)循环结构69【例8】 求100,1000以内的全部素数。 n=0; for m=100:1000 flag=1; j=m-1; i=2; while i solve(a*x2+b*x+c=0)u 求的根求的根 f (x) = (cos x)2 的一次导数的一次导数 x=sym(x); diff(cos(x)2)u 计算计算 f (x) = x2 在区间在区间

24、a, b 上的定积分上的定积分 syms a b x; int(x2,a,b)102q 在进行符号运算时，必须先定义基本的符号对象，可以是符号常量、符号变量、符号表达式等。符号对象是一种数据结构。符号对象与符号表达式q 含有符号对象的表达式称为符号表达式，Matlab 在内部把符号表达式表示成字符串，以与数字变量或运算相区别。q 符号矩阵/数组：元素为符号表达式的矩阵/数组。103u sym 函数用来建立单个符号变量，一般调用格式为：q 符号对象的建立：符号对象的建立：sym 和和 syms符号对象的建立例： a=sym(a) 符号变量 = sym(A)参数 A 可以是一个数或数值矩阵，也可以

25、是字符串a 是符号变量b 是符号常量 b=sym(1/3)C 是符号矩阵 C=sym(1 ab; c d)104q 符号对象的建立：符号对象的建立：sym 和和 syms符号对象的建立u syms 命令用来建立多个符号变量，一般调用格式为：syms 符号变量符号变量1 符号变量符号变量2 . 符号变量符号变量n 例： syms a b c a=sym(a); b=sym(b); c=sym(c);105q 符号表达式的建立：符号表达式的建立：例：建立符号表达式通常有以下2种方法：(1) 用 sym 函数直接建立符号表达式。(2) 使用已经定义的符号变量组成符号表达式。 y=sym(sin(x)

26、+cos(x) x=sym(x); y=sin(x)+cos(x)符号表达式的建立 syms x; y=sin(x)+cos(x)106Matlab 符号运算采用的运算符和基本函数，在形状、名称和使用上，都与数值计算中的运算符和基本函数完全相同符号对象的基本运算q 基本运算符基本运算符u 普通运算：+ - * / u 数组运算：.* . ./ .u 矩阵转置： .例： X=sym(x11,x12;x21,x22;x31,x32); Y=sym(y11,y12,y13;y21,y22,y23); Z1=X*Y; Z2=X.*Y;107符号对象的基本运算sin、cos、tan、cot、sec、cs

27、c、asin、acos、atan、acot、asec、acsc、exp、log、log2、log10、sqrtabs、conj、real、imagrank、det、inv、eig、lu、qr、svddiag、triu、tril、expm三角函数与反三角函数、指数函数、对数函数等q 基本函数基本函数108q 查找符号表达式中的符号变量查找符号表达式中的符号变量若表达式中有两个符号变量与 x 的距离相等，则ASCII 码大者优先。查找符号变量findsym(expr)按字母顺序列出符号表达式 expr 中的所有符号变量findsym(expr, N)按顺序列出 expr 中离 x 最近的 N 个符

28、号变量常量 pi, i, j 不作为符号变量109符号表达式的替换subs(f,x,a) 用用 a 替换字符函数替换字符函数 f 中的字符变量中的字符变量 x a 是可以是是可以是 数数/数值变量数值变量/表达式表达式 或或 字符变量字符变量/表达式表达式若 x 是一个由多个字符变量组成的数组或矩阵，则 a 应该具有与 x 相同的形状的数组或矩阵。q 用给定的用给定的数据数据替换符号表达式中的指定的替换符号表达式中的指定的符号变量符号变量110subs 举例 f=sym(2*u); subs(f,u,2) f2=subs(f,u,u+2) a=3; subs(f2,u,a+2) subs(f2

29、,u,a+2) syms x y f3=subs(f,u,x+y) subs(f3,x,y,1,2)ans=4f2=2*(u+2)ans=14ans=2*(a+2)+2)f3=2*x+2*yans=6u 例：指出下面各条语句的输出结例：指出下面各条语句的输出结果果f=2*u111符号矩阵符号矩阵 A=sym(1+x, sin(x); 5, exp(x)u 使用 sym 函数直接生成u 将数值矩阵转化成符号矩阵u 符号矩阵中元素的引用和修改 B=2/3, sqrt(2); 5.2, log(3); C=sym(B) A=sym(1+x, sin(x); 5, exp(x); A(1,2) % 引

30、用引用 A(2,2)=sym(cos(x) % 重新赋值重新赋值112六类常见符号运算六类常见符号运算q 因式分解、展开、合并、简化及通分等因式分解、展开、合并、简化及通分等q 计算极限计算极限q 计算导数计算导数q 计算积分计算积分q 符号求和符号求和q 代数方程和微分方程求解代数方程和微分方程求解113计算极限limit(f,x,a): 计算计算limit(f,a): 当当默认变量默认变量趋向于趋向于 a 时的极限时的极限limit(f): 计算计算 a=0 时的极限时的极限limit(f,x,a,right): 计算右极限计算右极限limit(f,x,a,left): 计算左极限计算左极

31、限lim( )xaf x例：计算 ，hxhxLh)ln()ln(lim0nnnxM1lim syms x h n; L=limit(log(x+h)-log(x)/h,h,0) M=limit(1-x/n)n,n,inf)114计算导数g=diff(f,v)：求符号表达式求符号表达式 f 关于关于 v 的导数的导数g=diff(f)：求符号表达式求符号表达式 f 关于关于默认变量默认变量的导数的导数g=diff(f,v,n)：求求 f 关于关于 v 的的 n 阶导数阶导数q diff syms x; f=sin(x)+3*x2; g=diff(f,x)115计算积分int(f,v,a,b):

32、计算定积分计算定积分int(f,a,b): 计算关于计算关于默认变量默认变量的定积分的定积分int(f,v): 计算不定积分计算不定积分int(f): 计算关于计算关于默认变量默认变量的不定积分的不定积分 syms x; f=(x2+1)/(x2-2*x+2)2; I=int(f,x) K=int(exp(-x2),x,0,inf)( )baf v dv ( )f v dv 例：计算 和2221(22)xIdxxx 20 xKedx 116符号求和 syms n; f=1/n2; S=symsum(f,n,1,inf) S100=symsum(f,n,1,100)symsum(f,v,a,b)

33、: 求和求和symsum(f,a,b): 关于关于默认变量默认变量求和求和( )bv af v 例：计算级数 及其前100项的部分和211nSn 例：计算函数级数21nxSn syms n x; f=x/n2; S=symsum(f,n,1,inf)117代数方程求解solve(f,v)：求方程关于指定自变量的解，求方程关于指定自变量的解，f 可以是可以是用字符串表示的方程用字符串表示的方程、符号表达式符号表达式或或符号方程符号方程；l solve 也可解方程组也可解方程组(包含非线性包含非线性)；l 得不到解析解时，给出数值解。得不到解析解时，给出数值解。求方程的根求方程的根从现代数学的眼光

34、来看从现代数学的眼光来看，傅里叶变换是一种特殊的积分傅里叶变换是一种特殊的积分变换。它能将满足一定条件的某个函数表示成正弦基函变换。它能将满足一定条件的某个函数表示成正弦基函数的线性组合或者积分。数的线性组合或者积分。傅立叶原理表明：任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加。而根据该原理创立的傅立叶变换算法利用直接测量到的原始信号，以累加方式来计算该信号中不同正弦波信号的频率、振幅和相位。0.506-0.506127“一只蝴蝶在巴西扇动翅膀，一只蝴蝶在巴西扇动翅膀，便有可能会在美国的德克萨斯便有可能会在美国的德克萨斯引起一场龙卷风引起一场龙卷风”。蝴蝶效应蝴蝶效应

35、E.N. E.N. 洛伦兹洛伦兹1972/12/29 1972/12/29 ，美国科学发展学，美国科学发展学会第会第 139 139 次会议次会议对非线性问题的求解，通常遇到常微分方程对非线性问题的求解，通常遇到常微分方程(组）或偏微分方程组）或偏微分方程(组），组），复杂微分方程很多情况下只有数值解。复杂微分方程很多情况下只有数值解。第五章 解常微分方程一、算法原理：龙格一、算法原理：龙格-库塔法库塔法以这种方法为基础，教材后续有精度要求不同的三阶、四阶龙格以这种方法为基础，教材后续有精度要求不同的三阶、四阶龙格库塔法原理介绍（课后阅读）库塔法原理介绍（课后阅读）第五章第五章 解常微分方程解

36、常微分方程二、解常微分方程的二、解常微分方程的Matlab指令及其用法指令及其用法解法的思想：二阶化为一阶，再用龙格解法的思想：二阶化为一阶，再用龙格-库塔法解一阶微分方程库塔法解一阶微分方程Matlab求解常微分方程基本步骤求解常微分方程基本步骤例题例题1：设非线性振子受 和 两个强迫力的作用，振子的运动方程为:tF11costF22costFtFxkxkdtxd221132122coscos试画出它的振动图像图，并讨论这种振动对初始值的影响？试画出它的振动图像图，并讨论这种振动对初始值的影响？解答：设：dtdxyxy则上述方程可以简化为21;tFtFykykdtdyydtdycoscos3

37、1211221假设：8 . 0, 1 . 0, 1, 5 . 0212121kkFF初始值：(1); 1 . 0|, 0000ttdtdxvx;1001. 0|, 0000ttdtdxvx (2)tFtFykykdtdyydtdycoscos31211221自定义函数zhendong.mfunction ydot=zhengdong(t,y,flag,k1,k2,F1,F2,w1,w2)ydot=y(2);-k2*y(1)3+k1*y(1)+F1*cos(w1*t)+F2*cos(w2*t);主函数：主函数：F1=0.5;F2=0.5;w1=0.5;w2=0.5;k1=1;k2=0.8;x0=

38、0.1;t1,x1=ode45(zhendong,0:0.01:100,x0,0.1,k2,k1,w1,w2,F1,F2);t2,x2=ode45(zhendong,0:0.01:100,x0,0.1001,k2,k1,w1,w2,F1,F2);plot(t1,x1,b,t2,x2,r)计算结果：蓝线表示初速度为1.001m/s,红线表示初速度为1.000m/s例题例题2-单摆的周期与摆角关系研究：大角度摆动与小角度摆动单摆的周期与摆角关系研究：大角度摆动与小角度摆动 单摆的摆动周期到底与摆角是否无关？单摆运动方程： 0,sin50sin2222lgdtdlgdtdo方程可以简化为时有摆角小于

39、假设:12121sin; 2;ylgdtdyydtdy：dtdyy振动方程可以写为clear alll=1;g=9.8;theta0=0.0873;t=(0:0.01:2)*pi;t1,u1=ode45(ex3313f,0:0.01:10,pi/7,0,);t2,u2=ode45(ex3313f,0:0.01:10,pi/3,0,);t3,u3=ode45(ex3313f,0:0.01:10,pi/10,0,);t4,u4=ode45(ex3313f,0:0.01:10,pi/20,0,);plot(t1,u1(:,1),r-.,t2,u2(:,1),t,0,t3,u3(:,1),m,t4,u

40、4(:,1),g),legend (pi/7,pi/3,x,pi/10,pi/20);hold onfunction dydt=ex3313f(t,y,flag);dydt=y(2);-9.8*sin(y(1);程序运行结果程序运行结果结论：在小角度摆动范围内，周期与振幅无关!Matlab中最小的数eps=2.2204e-016 5.10 5.10 边值问题和本征值问题边值问题和本征值问题一、问题描述bxayyxfy),( 第一类边值条件：第二类边值条件：第三类边值条件：.)(,)(byay.)(,)(byay.)()(,)()(1010bybyayay基本思想：将边值问题化为初值问题，在任一

41、边界上补充一个将边值问题化为初值问题，在任一边界上补充一个猜测的边界条件按初值问题求解方程，通常需要多次修正边界猜测的边界条件按初值问题求解方程，通常需要多次修正边界条件或猜测解，重新解方程，直到找出解为止。条件或猜测解，重新解方程，直到找出解为止。例：打靶法求弦振动方程222kdtd0) 1 ()0(m=3;global k;k=1;tol=1e-8;for n=1:m dk=k/15;k=k+dk; x,phi=ode45(ch5x2fun,0,1,0,1e-3); oldphi=phi(end,1); dphi=oldphi; while abs(dphi)tol k=k+dk; x,p

42、hi=ode45(ch5x2fun,0,1,0,1e-3); dphi=phi(end,1); if dphi*oldphi0 k=k-dk;dk=dk/2; end endplot(x,phi(:,1); hold on; kk(n)=k;endfunction yy=ch5x2fun(x,phi)global kyy=phi(2);-k2*phi(1);solinit=bvpinit(x,v,parameters)生成生成bvp4c调用指令所必须的调用指令所必须的“猜测解猜测解”sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,)给出微分方程边值问题

43、的近似解给出微分方程边值问题的近似解sxint=deval(sol，xint)计算微分方程积分区间内任何一点的解值计算微分方程积分区间内任何一点的解值二、用二、用Matlab指令求解边值问题指令求解边值问题solinit=bvpinit(x,v,parameters) x指定边界区间指定边界区间a,b上的初始网络，通常是等距排列的上的初始网络，通常是等距排列的（1M）一维数组。）一维数组。注意：使注意：使x(1)=a，x(end)=b；格点要单；格点要单调排列。调排列。 v是对解的初始猜测是对解的初始猜测 solinit（可以取别的任意名）带如下两个域：（可以取别的任意名）带如下两个域：sol

44、init.x是表示初始网格有序节点的（是表示初始网格有序节点的（1M）一维数组，）一维数组，并且并且solinit.x(1)一定是一定是a，solinit.x(end)一定是一定是b。M不不宜取得太大，宜取得太大，10数量级左右即可。数量级左右即可。solinit.y是表示网点上微分方程解的猜测值的（是表示网点上微分方程解的猜测值的（NM）二维数组。二维数组。solinit.y(:,i)表示节点表示节点solinit.x(i)处的解的猜处的解的猜测值。测值。初始解生成函数：初始解生成函数：bvpinit()sol=bvp4c(odefun,bcfun,solinit,options,p1,p2

45、,)输入参数：输入参数： odefun是计算导数的是计算导数的m函数文件。该函数的基本形式为：函数文件。该函数的基本形式为：dydx=odefun(x,y,parameters,p1,p2,)，在此，自变量，在此，自变量x是标量，是标量，y，dydx是列向量。是列向量。 bcfun是计算边界条件的是计算边界条件的m函数文件。其基本形式为：函数文件。其基本形式为：res=bcfun(ya,yb,parameters,p1,p2,)， 输入宗量输入宗量options是用来改变是用来改变bvp4c算法的控制参数的。算法的控制参数的。在最基本用法中，它可以缺省，此时一般可以获得比较满在最基本用法中，它

46、可以缺省，此时一般可以获得比较满意的边值问题解。如需更改可采用意的边值问题解。如需更改可采用bvpset函数，使用方法函数，使用方法同同odeset函数。函数。 输入量输入量p1,p2等表示希望向被解微分方程传递的已知参数。等表示希望向被解微分方程传递的已知参数。如果无须向微分方程传递参数，它们可以缺省。如果无须向微分方程传递参数，它们可以缺省。输出参数：输出参数：输出变量输出变量sol是一个结构体是一个结构体usol.x是指令是指令bvp4c所采用的网格节点；所采用的网格节点；usol.y是是y(x)在在sol.x网点上的近似解值；网点上的近似解值；usol.yp是是y(x)在在sol.x网

47、点上的近似解值；网点上的近似解值；usol.parameters是微分方程所包含的未知参数的近似解值。是微分方程所包含的未知参数的近似解值。 当被解微分方程包含未知参数时，该域存在。当被解微分方程包含未知参数时，该域存在。边值问题求解指令：边值问题求解指令：bvp4c()求解边值问题举例求解边值问题举例 例：求二阶方程 满足边界条件 的解。 0 zcz2)4(, 0)0(zz（1）首先把原方程改写为一阶方程组，原边界条件改写为“残数”形式：令y1=z, y2=z, 方程变化为同时，边界条件变为：ya(1)=0, yb(1)=-2. 即：bc(ya(1),yb(1)+2)=0(2)(1)(2)(

48、1)yyyycy （2）根据改写后的方程和边界条件编写相应的函数M文件：twoode.mfunction dydx=twoode(x,y,c)dydx= y(2) , -c*abs(y(1); %导数列向量twobc.mfunction res = twobc(ya,yb,c) % ya,yb分别为左右边界处的解变量；位置不可改变，变量名可以任取 res= ya(1), yb(1)+2; %边界残差列向量（3）猜测解创建：令初始网格为linspace(0,4,5),初始猜测解向量为1;0,于是猜测解由指令sinit求得：sinit=bvpinit(linspace(0,4,5),1;0); %

49、5个均匀分布网格节点上个均匀分布网格节点上赋初始解赋初始解（4）求指定网格上的解：c=1;sol=bvp4c(twoode,twobc,sinit, ,c);%属性选项为属性选项为，取缺省，取缺省设置设置 （5）求积分区间内任意点上的解：x=linspace(0,4,100);y = deval(sol,x); % bvpval(sol,x); plot(x,y(1,:),b-,sol.x,sol.y(1,:),ro,sinit.x,sinit.y(1,:),ks)legend(fontname隶书隶书fontsize20插值后的解曲线插值后的解曲线,解点解点,猜猜测解点测解点,0) 说明：t

50、woode.m中的导数变量和twobc.m中的残差向量一定要写成列向量形式。原方程组等价于以下标准形式的原方程组等价于以下标准形式的方程组：方程组：solinit=bvpinit(linspace(0,1,10),1 0);sol=bvp4c(ODEfun,BCfun,solinit);x=0:0.05:0.5;y=deval(sol,x);xP=0:0.1:0.5;yP=0 -0.41286057 -0.729740656. -0.95385538 -1.08743325 -1.13181116;plot(xP,yP,o,x,y(1,:),r-)legend(Analytical Solut

51、ion,Numerical Sol.)% 定义定义ODEfun函数函数function dydx=ODEfun(x,y)dydx=y(2);y(1)+10;% 定义定义BCfun函数函数function bc=BCfun(ya,yb)bc=ya(1);yb(1);10(0)(1)0yyyy求解两点边值问题：求解两点边值问题： 121,yy yy令：令： 122110 yyyy11(0)0,(1)0yy边界条件为：边界条件为： 边值问题的求解举例边值问题的求解举例原方程组等价于以下标准形式的原方程组等价于以下标准形式的方程组：方程组：solinit=bvpinit(linspace(0,1,10

52、),0 1);sol=bvp4c(ODEfun,BCfun,solinit);x=0:0.1:1;y=deval(sol,x);plot(x,y(1,:),r-)legend(Analytical Solution,Numerical Solution,.location,Northwest)legend boxoff% 定义定义ODEfun函数函数function dydx=ODEfun(x,y)dydx=y(2);(1+x2)*y(1)+1;% 定义定义BCfun函数函数function bc=BCfun(ya,yb)bc=ya(1)-1;yb(1)-3;求解：求解： 121,yy yy令

53、：令： 边界条件为：边界条件为： 2(1)1(0)1, (1)3yxyyy12221(1)1 yyyxy03) 1 (01)0(11yy边值问题的求解举例边值问题的求解举例c=1;solinit=bvpinit(linspace(0,4,10),1 1);sol=bvp4c(ODEfun,BCfun,solinit,c);x=0:0.1:4;y=deval(sol,x);plot(x,y(1,:),b-,sol.x,sol.y(1,:),ro)legend(解曲线解曲线,初始网格点解初始网格点解)% 定义定义ODEfun函数函数function dydx=ODEfun(x,y,c)dydx=y

54、(2);-c*abs(y(1);% 定义定义BCfun函数函数function bc=BCfun(ya,yb,c)bc=ya(1);yb(1)+2;求解：求解： 121,zz zz令：令： 0 zc z2)4(, 0)0(zz边值问题的求解举例边值问题的求解举例1csolinit=bvpinit(linspace(0,pi,10),1;1,lmb);opts=bvpset(Stats,on);sol=bvp4c(ODEfun,BCfun,solinit,opts);lambda=sol.parametersx=0:pi/60:pi;y=deval(sol,x);plot(x,y(1,:),b-,sol.x,sol.y(1,:),ro)legend(解曲线解曲线,初始网格点解初始网格点解)% 定义定义ODEfun函数函数function dyd