数学物理建模作业讲评.ppt

《数学物理建模作业讲评.ppt》由会员分享，可在线阅读，更多相关《数学物理建模作业讲评.ppt（38页珍藏版）》请在三一办公上搜索。

1、对函数sin(x)+cos(x)取任意多个数据点，求拟合6阶多项式，并图示拟合情况，并对图形进行修饰。,x=0:0.1:20;y=sin(x)+cos(x);p=polyfit(x,y,6);f=polyval(p,x);plot(x,y,o,x,f,-);,用符号运算应该怎么做？%计算位移和速度function s,v=velocity(X0,V0,A,T);syms x x0 v0 a ty1=dsolve(D2x=a,Dx(0)=v0,x(0)=x0,t);%求解满足的微分方程，得到位移y2=diff(y1,t);%对位移作一次导数得到速度s=subs(y1,x0,v0,a,t,X0,V

2、0,A,T);%变量替换得到位移和速度的值s=double(s);v=subs(y2,x0,v0,a,t,X0,V0,A,T);v=double(v);,求点P(2,-4,5)关于直线的垂足和对称点。要求编M函数文件来实现，并绘制这些直线和点。,%解法1：求p点垂足和对称点syms x y zp=2,-4,5;q=x,y,z;g1=(x-2)/3=(y+3)/7;g2=(y+3)/7=(z-1)/(-3);g3=dot(q-p,3 7-3);f=solve(g1,g2,g3)C=double(f.x,f.y,f.z)S=2.*C-p,垂足点(1.1493,-4.9851,1.8507),对称点

3、(0.2985,-5.9701,-1.2985),%求中点和对称点%入口参数 P0 已知点坐标%x0 直线上点坐标%k 直线方向矢量function C,S=points(P0,x0,k);%x0=2,-3,1;%k=3 7-3;%P0=2,-4,5;P1=sym(x,y,z);%符号表达式，待求对称点坐标P2=(P0+P1)/2;%符号表达式，垂足点坐标Online=(P2-x0)./k;%符号表达式，垂足点代入直线方程Perpen=sum(P2-x0).*k);%符号表达式，垂直条件代入f1=sym(strcat(char(Online(1),=,char(Online(2);%垂足点在线

4、条件1f2=sym(strcat(char(Online(1),=,char(Online(3);%垂足点在线条件2f3=sym(strcat(char(Online(2),=,char(Online(3);%垂足点在线条件3f4=sym(strcat(char(Perpen),=0);%连线垂直于直线条件s=solve(f1,f2,f4);%求条件方程，注意这里4个条件只有3个是独立的，求解3个未知数S=double(s.x,s.y,s.z);C=double(S+P0)/2);,解法2：function Pd,Pp=fp(P,Vp,V)Pv=P-Vp;L=sqrt(dot(V,V),C=V

5、*dot(Pv,V)/L*2D=Pv-C;Pd=Vp+C;Pp=P-2*D;Pd,Pp=fp(2,-4,5,2,-3,1,3,7,-3),解法3：k=3,7,-3;p=2,-4,5;A=7-3 0;0 3 7;3 7-3;b=23-2-37;x=inv(A)*bm=2*x-p,回顾一下泰勒级数展开泰勒(Taylor)展开定理设f(z)在区域D:|z-z0|R内解析，则在D内f(z)可展开成泰勒级数其中,一些常见的泰勒级数展开公式,回顾一下罗朗级数展开罗朗级数展开定理设函数f(z)在D:R1|z-z0|R2内解析，则在此圆环内f(z)必可展开成罗朗级数其中,应用两种级数展开复变函数的时候首先确定

6、：展开所在的区域属于哪一个环域若内部区域解析，则采用泰勒展开若内部有奇点或不解析，则采用罗朗展开,绘制复变函数 的图像，并验证其泰勒展开和罗朗展开,泰勒展开,罗朗展开,罗朗展开,图形的绘制(1)直接绘制z=4*cplxgrid(40);cplxmap(z,1./(z+eps*(abs(z)=1)-1)./(z+eps*(abs(z)=2)-2),10*pi);,图形的绘制(2)泰勒展开 z=2*cplxgrid(30);z1=z;z1(abs(z1)=1)=NaN;w1=1;u1=1;p1=1;q1=1;for k=1:100 u1=u1.*z1;w1=w1+u1;p1=p1.*z1/2;q1

7、=q1+p1;endcplxmap(z1,w1+q1/2)colorbar,图形的绘制(3)罗朗展开z=3*cplxgrid(30);z2=z;z2(abs(z2)=2)=NaN;w2=1./z2;u2=1./z2;p2=1;q2=1;for k=1:100 u2=u2./z2;w2=w2+u2;p2=p2.*z2/2;q2=q2+p2;endcplxmap(z2,-w2-q2/2)colorbar,图形的绘制(4)罗朗展开z=3*cplxgrid(30);z2=z;z2(abs(z2)=2)=NaN;w2=1./z2;u2=1./z2;p2=2./z2;q2=2./z2;for k=1:10

8、0 u2=u2./z2;w2=w2+u2;p2=2*p2./z2;q2=q2+p2;endcplxmap(z2,-w2+q2/2)colorbar,直接用Matlab怎么绘制图形？z=2*cplxgrid(30);w=1./(z+eps*(z=0).2./(z-i+eps*(z=i);w=1./(z.2+eps*(z=0)./(z-i+eps*(z=i);w=1./(z.2.*(z-i)+eps*(z=0)+eps*(z=i);cplxmap(z,w,10*pi);colorbar(vert);,n=input(input n:=);%取展开式前几项项数syms t z z1c=1:n;d=1

9、:n;for k=1:n z=(1/2)*cos(t)+i*(1/2)*sin(t)+i;%在0|z-i|1区域内展开，令z-i=R*cos(t)+i*R*sin(t),0R1,可取R=1/2 f=1/(z2)*(z-i);%函数表达式 c(k)=double(1/(2*pi*i)*int(f/(z-i)(k+1)*diff(z),t,0,2*pi);%求罗朗级数正幂项系数endfor m=-n:-1 z=(1/2)*cos(t)+i*(1/2)*sin(t)+i;f=1/(z2)*(z-i);d(m+n+1)=double(1/(2*pi*i)*int(f/(z-i)(m+1)*diff(z

10、),t,0,2*pi);%求罗朗级数负幂项系数endq=double(1/(2*pi*i)*int(f/(z-i)(0+1)*diff(z),t,0,2*pi);%求罗朗级数0次项系数input(f的罗朗级数的(z-i)-n项到(z-i)n项的系数依次是：)d,q,c,运行结果：ans=Columns 1 through 6 0 0 0 0-1.0 0-2.0i Columns 7 through 11 3.0 0+4.0i-5.0 0-6.0i 7.0000,z=2*cplxgrid(30);z1=z-i;z1(abs(z1)=1|abs(z1)=0)=NaN;w1=0;u1=i.*(z1+

11、eps).(-2);for k=-1:10 u1=u1.*i.*(z1+eps);w1=w1+(k+2).*u1;endcplxmap(z,w1,10*pi)colorbar(vert),z=2*cplxgrid(30);z2=z-i;z2(abs(z1)=1)=NaN;w2=0;u2=i.*(z2+eps).(-2);for k=3:13 u2=u2./i./(z2+eps);w2=w2+(k-2).*u2;endcplxmap(z,w2,10*pi)colorbar(vert),w1(find(isnan(w1)=1)=0;w2(find(isnan(w2)=1)=0;w=w1+w2;cp

12、lxmap(z,w,10*pi)colorbar(vert),%cplxgrid.mfunction z=cplxgrid(m)r=(0:m)/m;theta=pi*(-m:m)/m;z=r*exp(i*theta);,%cplxroot.mfunction cplxroot(n,m)r=(0:m)/m;theta=pi*(-n*m:n*m)/m;z=r*exp(i*theta);s=r.(1/n)*exp(i*theta/n);surf(real(z),imag(z),real(s),imag(s);,绘制极坐标系下的网格z=cplxgrid(30);theta=angle(z);r=abs

13、(z);polar(theta,r,k),另一种做法m=20;r=linspace(0,1,m+1);theta=linspace(-pi,pi,2*m+1);R,Q=meshgrid(r,theta);X,Y=pol2cart(Q,R);Z=ones(size(X);mesh(X,Y,Z)view(0,90)axis equalaxis off,验证Laplace时移性质,syms t st0=sym(t0,positive)syms t0 positiveft=sym(f(t-t0)*sym(Heaviside(t-t0)FS=laplace(ft,t,s),IFS=ilaplace(FS

14、,s,t),FS：exp(-s*t0)*laplace(f(t),t,s)IFS：f(t-t0)*heaviside(t-t0),求Fourier变换，x是参量，t是变量。,syms x t vf=exp(-(t-x)*heaviside(t-x)g=fourier(f,t,v),g：1/(1+i*v)*exp(-i*x*v),V(x,y)=log(x2+y2);V=log(x.2+y.2+eps);xMax=10;yMax=10;NGrid=20;xPlot=linspace(-xMax,xMax,NGrid);x,y=meshgrid(xPlot);VPlot=eval(V);ExPlot,EyPlot=-gradient(VPlot);gradient(-VPlot)subplot(1,2,1),meshc(VPlot);xlabel(x);ylabel(y);zlabel(电位)subplot(1,2,2)axis(-xMax,xMax,-yMax,yMax);cs=contour(x,y,VPlot);hold onquiver(x,y,ExPlot,EyPlot);xlabel(x);ylabel(y);,?Error using=uminusToo many output arguments.,