课堂授课专题1：数学物理建模思想与数据可视化方法.ppt

数学物理建模与计算机辅助设计,专题1：数学物理建模思想与 数据可视化方法,江萍 联系电话QQ：397379519 邮箱： 课程QQ群： 245625386,Page 2,上课时间及地点,Page 3,Page 4,课程基本信息,先修课程：数学物理方法 课程教材： 数学物理方程的Matlab解法与可视化 彭芳麟，清华大学出版社 参考教材： 杨华军 数学物理方法与仿真 薛定宇 高等应用数学问题的MATLAB求解 主要内容: 留数和级数的计算 复变函数的积分、变换和逆变换 特殊函数的可视化 数学物理方程的求解,课程安排,Page 5,专题1：数学物理建模思想和数据可视化,专题2：复变函数的计算机仿真,专题3：特殊函数的可视化,专题4：数学物理方程的求解与可视化,专题5：数学物理方程-稳定场方程的求解与可视化 数学物理方程-热传导方程的求解与可视化,专题6：空间光通信实践与数学物理建模,专题4：数学物理方程-波动方程的求解与可视化,专题3：特殊函数的可视化,专题1：MATLAB基础和数据可视化,专题2：复变函数的计算机仿真,理论授课,（每个专题2学时）,（每个专题4学时）,上机实验,Page 6,课程的重点难点,重点： 计算机仿真绘图 数学物理方程的求解与可视化 难点： 数学物理方程的求解与可视化 特点： 计算机编程能力和程序设计能力,Page 7,课程考核,考核形式：考查 考核内容和比重： 完成课后作业 20% 考勤情况 20% 课程设计（大作业） 60% 特别提示： 严禁抄袭，否则相应成绩(课后作业或大作业)记0 严禁上机时上网、玩游戏、聊QQ等，否则上机成绩记0,Page 8,课程授课计划和安排,提交作业邮箱： 397379519 作业提交形式： 文件格式：WORD文件(*.doc或*.docx) 文件命名规则：学号+姓名.doc(或.docx) 如：第一次作业_201404021001 _张三. doc 或：第一次作业_ 201404021001 _张三. docx,（一）数学物理方法问题建模的基本原理,Page 9,主要解决的问题： （1）三维可视化方法使得“数学物理方法”中复杂、抽象、烦琐的理论模型变得直观、具体、明了。,例1: 绘制复变函数 的图像，并验证其泰勒展开和罗朗展开的函数。,图1 收敛圆解析示意图,图2 在|z|1内泰勒级数展开图形,（一）数学物理方法问题建模的基本原理,Page 10,图3 在1|z|2内罗朗级数展开图形,图4 |z|2内罗朗级数展开图形,（一）数学物理方法问题建模的基本原理,Page 11,图5 特殊函数复变量(z)函数的图形,例2: 绘制复变量(z)函数的图形,（一）数学物理方法问题建模的基本原理,图6 球函数 的图形,Page 12,例3： 绘制球函数 的图形，取l=3，m=0,（一）数学物理方法问题建模的基本原理,Page 13,图7 矩形膜振动的前四个本征函数三维可视化结果,例4：绘制四周固定的矩形膜本征振动的图形,（一）数学物理方法问题建模的基本原理,Page 14,（2）通过仿真设计与验证，解决前沿性的科研中缺乏实验条件的问题。,图9 布拉格光子晶体光纤场振幅Hz三维分布,图8 布拉格光纤横截面结构示意图,主要解决的问题：,Page 15,Page 16,Page 17,（二）利用MATLAB进行数学物理建模基础,例1：已知函数，如何求导及高阶导数？,问题：求导过程很繁杂，容易出错,怎么办?,思路： 由分式求导公式，得出, 逐次求导则可以得出,使用Matlab的符号运算功能 syms x diff(sin(x)/(x2+4*x+3),4) ans= sin(x)/(x2+4*x+3)+4*cos(x)/(x2+4*x+3)2*(2*x+4)-12*sin(x)/(x2 +4*x+3)3*(2*x+4)2+12*sin(x)/(x2+4*x+3)2-24*cos(x)/(x2+ 4*x+3)4*(2*x+4)3+48*cos(x)/(x2+4*x+3)3*(2*x+4)+24*sin(x)/(x2+4*x+3)5*(2*x+4)4-72*sin(x)/(x2+4*x +3)4*(2*x+4)2+ 24*sin(x)/(x2+4*x+3)3,Page 18,MATLAB的常量和变量,Matlab预定义一些常量,Page 19,Page 20,MATLAB的向量的生成,向量的生成的三种方式 1 直接输入向量 a=1 2 3 4 5 6 7 8 9 b=1;2;3;4;5;6;7;8;9 用空格或逗号生成行向量,用分号生成列向量 2 利用冒号表达式基本形式:x=x1:step:xn a=1: 2: 12 a=1 3 5 7 9 11 3 线性等分向量生成 基本格式y=linspace(x1,x2,n) y=linspace(0,100,6) y=0 20 40 60 80 100,Page 21,MATLAB的向量的运算,加减与数加减运算 A+B A-B A+b A-b 加减规则：(1)长度相同；(2)对应元素加减 数加减规则：每个元素都加减同一个数 数乘除运算 A*b A/b 运算规则：每个元素都乘除以同一个数 带点乘除运算 A.*B A./B 运算规则： (1)长度相同；(2)对应元素相乘除 点积和叉积运算 点积： c=dot(a,b) 返回A和B的数量点积，A和B须同维度 叉积： c=cross(a,b) 返回A和B叉积，A和B维度必须同为3,A=1,2,3;B=3,4,5; C=dot(A,B); D=cross(A,B);,Page 22,MATLAB的矩阵的生成,生成矩阵的四种方式 直接输入小矩阵 A=1 2 3; 4 5 6; 7 8 9 利用特殊矩阵生成函数生成矩阵 利用已有矩阵构建新矩阵 B=A(1:2,:) C=repmat(B,2,3),Page 23,MATLAB的矩阵的运算,矩阵的四则运算 矩阵加减运算 + - 规则：矩阵同阶(维数相同，各维度长度对应相同)，对应做加减 A=1 2 3;2 3 4;3 4 5; B=1 1 1;2 2 2;3 3 3; C=A+B;,Page 24,MATLAB的矩阵的运算,矩阵的四则运算 矩阵的乘法运算 * 规则：(1) 若A为ij阶，B必须为jk阶时A和B才能相乘 (2) A=1 2 3;2 3 4;3 4 5; B=1 1 1;2 2 2;3 3 3; E=B,5 5 5;C=A*E;,Page 25,MATLAB的矩阵的运算,矩阵的四则运算 矩阵的除法运算，分为左除和右除/运算 规则：(1)左除“: X=AB是方程式AX=B的解 (2)右除“/“: X=B/A是方程式XA=B的解 通常用矩阵的除法求解方程组的解 A=5 12 8;6 5 8;9 6 10; B=7;11;7; X=AB;,Page 26,MATLAB的矩阵的运算,矩阵的行列式运算 det(A) 规则：矩阵为方阵 C=det(A);,Page 27,MATLAB的矩阵的运算,矩阵的幂运算 An 和 A.n 规则：(1) An中A必须为方阵，相当于n个方阵A连乘 (2) A.n中A不必为方阵，结果为A中各元素取n次幂 D=1,2,3;4,5,6;7,8,9; E=D3; F=D.3;,Page 28,MATLAB的多项式运算,多项式的表示 p(x)=x3-5x2+6x-33 p=1 -5 6 -33 p=1 -5 6 -33; A=poly2sym(p) A = x3-5*x2+6*x-33 求多项式的根 roots roots(p) ans = 5.0939 -0.0470 + 2.5448i -0.0470 - 2.5448i 求多项式的乘除 conv(乘法) deconv(除法) p= 1 2 3; q=4 5 6 7; pq=conv(p,q); deconv(pq,q); pq=4 13 28 34 32 21 ans=1 2 3 求多项式的微分 polyder Dpq= polyder(pq) Dpq=20 52 84 68 32,Page 29,MATLAB的符号计算功能,符号运算与数值运算的区别 数值运算必须先对变量赋值，再运算 符号运算运算结果以标准的符号形式（公式）表达 符号运算的特点 运算对象可以是没赋值的符号变量 可以获得任意精度的解,%浮点算数运算 1/2+1/3 ans = 0.83333333333333,%符号运算 sym(1/2+1/3) ans = 5/6 vpa(1/2+1/3,20) ans = .83333333333333333333,Page 30,MATLAB的符号计算功能,符号变量与符号表达式 符号变量 定义：syms a b c 规则：用空格隔开 符号表达式 定义： 字符串 或 符号方程,f1=a*x2+b*x+c 二次三项式 f2= a*x2+b*x+c=0 方程 f3=Dy+y2=1 微分方程,Page 31,MATLAB的符号计算功能,符号矩阵的创建：用Matlab函数sym创建矩阵 命令格式：A=sym( ) A = sym(a , 2*b ; 3*a , 0) A = a, 2*b 3*a, 0 符号矩阵的修改和转换 符号矩阵的修改：用A1=subs(A, new, old)来修改,例：A2=subs(A, c, b) A2 = a, 2*c 3*a, 0,Page 32,MATLAB的符号计算功能,符号矩阵与数值矩阵的转换 将数值矩阵转化为符号矩阵，函数调用格式：sym(A) 将符号矩阵转化为数值矩阵，函数调用格式：double(A),sym(A) ans = 1/3, 5/2 10/7, 2/5,double(A) ans = 0.3333 2.5000 1.4286 0.4000,A=1/3,2.5;1/0.7,2/5 A = 0.3333 2.5000 1.4286 0.4000,Page 33,MATLAB的符号计算功能,符号矩阵的运算 7.0以后版本对符号表达式直接用+ - * / 进行相应运算, f=sym(a,b;c,d) f = a, b c, d,g=sym(2*a,c;b,3*d) g = 2*a, c b, 3*d,f*g ans = 2*a2+b2, a*c+3*b*d 2*a*c+b*d, c2+3*d2,f.*g ans = 2*a2, c*b c*b, 3*d2,Page 34,MATLAB的符号计算功能,符号运算函数 inv 逆矩阵 transpose 矩阵的转置 simplify 符号简化运算 simple 符号矩阵简化, f=sym(cos(x)+sin(x); g=sym(cos(x)-sin(x); h=f*g h = (cos(x)+sin(x)*(cos(x)-sin(x) simplify(h) ans = 2*cos(x)2-1 simple(h) ans = cos(2*x),Page 35,MATLAB的符号计算功能,符号微积分 diff(f) 对缺省变量x求微分 diff(f,v) 对指定变量v求微分 diff(f,v,n) 对指定变量v求n阶微分 int(f) 对f表达式的缺省变量x求积分 int(f,v) 对f表达式的v变量求积分 int(f,v,a,b) 对f表达式的v变量在(a,b) 区间求定积分,F=int(int(x*exp(-x*y),x),y) F = 1/y*exp(-x*y),F=diff(sym(x3+3*x2*y+5*x*y2), x,2) F = 6*x+6*y,Page 36,MATLAB的符号计算功能,符号代数方程求解 solve 解一般的线性方程、非线性方程、代数方程和代数方程组 当方程组不存在符号解且无其他自由参数，则给出数值解,例1：f = ax2+bx+c 求解 f=a*x2+b*x+c; solve(f) %对缺省变量x求解 ans = 1/2/a*(-b+(b2-4*a*c)(1/2) 1/2/a*(-b-(b2-4*a*c)(1/2) solve(f , b) %对指定变量b求解 ans = -(a*x2+c)/x,Page 37,MATLAB的符号计算功能,例2：符号方程tan(2*x)=sin(x)求解 f1=solve(tan(2*x)=sin(x) f1 = acos(1/2 - 3(1/2)/2) acos(3(1/2)/2 + 1/2) 0 -acos(1/2 - 3(1/2)/2) -acos(3(1/2)/2 + 1/2) double(f1) ans = 1.9455 0 + 0.8314i 0 -1.9455 0 - 0.8314i,Page 38,MATLAB的符号计算功能,例3：解方程组x+y+z=1;x-y+z=2;2x-y-z=1; g1=x+y+z=1; g2=x-y+z=2; g3=2*x-y-z=1; f=solve(g1,g2,g3) f = x: 1x1 sym y: 1x1 sym z: 1x1 sym f.x, f.y, f.z ans= 2/3, -1/2, 5/6,Page 39,MATLAB的符号计算功能,符号微分方程求解 dsolve(f, g) f微分方程，可多至12个微分方程的求解；g为初始条件 微分方程的各阶导数项以大写字母D表示 Dy (y的一阶导数) D2y(y的二阶导数) y1,y2=dsolve(x1,x2,xn) 返回微分方程的解,例1：一阶微分方程 dsolve(Dx=y,Dy=-x,x(0)=0,y(0)=1) ans = x(t) = sin(t), y(t) = cos(t),例2：二阶 微分方程 dsolve(D2y=-a2*y,y(0)=1,Dy(pi/a)=0) ans = cos(a*x),Page 40,MATLAB的符号计算功能,例3：求下面方程的解,y=dsolve(D2y+2*Dy+2*y=0,y(0)=1,Dy(0)=0) y = exp(-x)*cos(x)+exp(-x)*sin(x) ezplot(y) %方程解y(t)的时间曲线图,Page 41,（三）MATLAB绘图和数据可视化,为何要将数据可视化 难以从大量数据或符号中感受具体含义 直观的图形更容易加深理解 快速、有效地表达想法、观察或结论 Matlab的数据可视化功能 点、线、面处理 二维、三维和四维表现图 图形着色、消隐、光照、渲染及多视角处理等,Page 42,Matlab的强大表现能力,Page 43,MATLAB基本图形绘制,二维直角坐标系图形绘制 plot(x,y) 以x为横坐标，y为纵坐标,x=0:0.01*pi:2*pi; y=sin(x); plot(x,y),Page 44,MATLAB基本图形绘制,双纵坐标系： plotyy,x=0:0.01*pi:2*pi; y=sin(x); z=exp(x); plotyy(x, y, x, z, plot, semilogy),绘图方式可以是 plot、semilogx、semilogy和loglog等,Page 45,MATLAB基本图形绘制,特殊坐标系的二维图形绘图 对数坐标系：semilogx semilogy loglog 极坐标系： polar,x=0.0:0.01*pi:4*pi; y=sin(x/2)+x; polar(x, y, -),Page 46,MATLAB基本图形绘制,特殊二维图形的绘图,1. 条形图 bar x=1:11; y=exp(-(x-6).2/4); bar(x, y),2. 离散点图 stem x = 0:0.25:(3*pi); stem(x,sin(x), sr, fill);,Page 47,MATLAB基本图形绘制,特殊二维图形的绘图,3. 饼图 pie x=2, 4, 6, 8; pie(x,math, English, Chinese, music),4. 等高线图 contour x=-pi:0.1*pi:pi; X,Y=meshgrid(x); Z=sin(X).*cos(Y); contour(Z),Page 48,MATLAB基本图形绘制,三维图形绘制 三维曲线绘制 plot3,例1：绘制三维螺旋线 x=0:pi/50:10*pi; y=sin(x); z=cos(x); plot3(x,y,z),例2：渐变的三维螺旋线 t=0:0.1:10*pi; x=exp(-t/20).*cos(t); y=exp(-t/20).*sin(t); z=t; plot3(x,y,z); xlabel(x); ylabel(y); zlabel(z);,Page 49,MATLAB基本图形绘制,三维图形绘制 三维曲面绘制的数据准备：创建3D网格和表面 meshgrid,x = -1 0 1; y = 9 10 11 12; X,Y = meshgrid(x,y) X =-1 0 1 -1 0 1 -1 0 1 -1 0 1 Y = 9 9 9 10 10 10 11 11 11 12 12 12,meshgrid(x,y)生成X，Y矩阵 X重复x为行向量 Y重复y为列向量,X,Y = meshgrid(linspace(0,2*pi,50),linspace(0,pi,50); Z = sin(X).*cos(Y); mesh(X,Y,Z) xlabel(x); ylabel(y); zlabel(z); axis(0 2*pi 0 pi -1 1),Page 50,MATLAB基本图形绘制,三维图形绘制 绘制三维曲面的函数 mesh 绘制网格曲面 surf 绘制完整曲面 meshc, surfc 带有轮廓线的表面图形 meshz 带有遮帘线（参考平面）的表面图形 surfl 源于指定方向的表面图形,t=0:pi/6:4*pi; x,y,z=cylinder(4+cos(t),30);