非线性动力系统-MatLab编程简介

MatLab编程简介MatLab的从例子都应该能看明白、先介绍MatLab里几个挺好用，但是大家以前可能没用到的命令，用过用法，所以我基本不加注释。1、求解初等方程：，+■=0solve('xA2+wA2')ans=i*w-i*w2、求解微分方程：X--2x,以及求函数的导数。xt=dsolve('D2x=-wA2*x','t'),yt=diff(xt,'t')xt=C1*sin(w*t)+C2*cos(w*t)yt=C1*cos(w*t)*w-C2*sin(w*t)*w3、符号作图：y二Asin(■t■)symst;%这个命令用来声明下面要用到的符号变量ezplot(sin(t+pi/6),[0,2*pi]);axis([0,2*pi,-1,1]);sin(t+1/6")t121224、符号作带等高线曲面图：zy22x222symsxyw;ezsurfc(y*y/2+2*x*x,[-4,4,-8,8]);5、画向量场P(x,y)二y,Q(x,y)二-sinxu=-4:0.3:4;[x,y]=meshgrid(u);%quiver(x,y,y,-sin(x));这个命令用来先生成一个数据网格以画向量场。6、极坐标作图：Theta=0:0.01:5*pi;_1_1-e^"polar(Theta,1./sqrt(1+exp(-2*Theta)),'r');7、画等高线:z二x2(xT)2y2symsxy;ezcontour(x*x*(x-1)*(x-1)+y*y,[-0.3,1.3,-0.5,0.5]);MatLab中查在线帮助。MatLab的一些目前可能用到的命令列表，需要知道更详细的格式，可以在ezmeshezmeshc画网格图画带等咼线的网格图ezsurf画曲面图ezsurfc画带等咼线的曲面图ezplot符号曲线图ezplot3三维符号曲线图ezpolar极坐标曲线图ezcontour画等咼线图ezcontourf画填充等咼线图plot数值二维图poar数值极坐标图plot3数值三维图mesh数值网格图meshc数值带等高线图surf数值曲面图surfc数值带等高线图quiver矢量图solve符号解代数方程dsolve符号解微分方程meshgrid产生数据网格gradient求数值梯度三、微分方程的数值解法微分方程u=f(t,u)，其中U为矢量。以下以简谐振动方程为例:y-_4x记Um=U(tm)1欧拉法：Umi=Umf(tm,Um)hh=0.001;N=2A16;t=O:h:N*h;x=zeros(size(t));y=zeros(size(t));%可以直接取x=zeros(1,N+1);但是一来容易范少算一个点的错误，二来下次改%度时，还得记得改x的长度，这样编程x的长度随时间变量t自动取比较好x(i)=i;y(i)=0;fori=1:N%这里的循环次数极容易弄错，要么多算要么少算一个点，这种错误往往还查不出来。x(i+1)=x(i)+y(i)*h;y(i+1)=y(i)-4*x(i)*h;end;Plot(x,y)t的长(1)2.521.510.50-0.5-1-1.5-2-0.5-2.5L-1.52、龙格库塔法：um1=um■匚hC2-2'2'4)6叫=f(tm,um)11.22羁=f(tm+丄h,Um+1⑷2"h)22$4=f(tm*h,Um+灼「h)h=0.001;h2=h/2;h6=h/6;N=2A19;%此处设置两个变量h2和h6是为了减少下面的大循环中的除法的次数，记住一个原则:%与循环变量无关，那么所有关于这个量的运算能算好的都先算好。t=0:h:N*h;x=zeros(size(t));y=zeros(size(t));x(1)=1;y(1)=0;fori=1:Nux=x(i);uy=y(i);wx1=uy;wy1=-4*ux;ux=x(i)+wx1*h2;uy=y(i)+wy1*h2;wx2=uy;wy2=-4*ux;(B.2)(B.3)如果在循环里有某个量ux=x(i)+wx2*h2;uy=y(i)+wy2*h2;wx3=uy;wy3=-4*ux;ux=x(i)+wx3*h;uy=y(i)+wy3*h;wx4=uy;wy4=-4*ux;x(i+1)=x(i)+(wx1+2*wx2+2*wx3+wx4)*h6;y(i+1)=y(i)+(wy1+2*wy2+2*wy3+wy4)*h6;end;%在上面的循环里，我设了两个临时变量ux,uy，目前大家能看到的是两个优点，一个是程序结构非常清楚,%另一个优点是这个程序的可移植性非常强，如果算另一个方程，只要作极小的改动就行了，见下一个例子%至于用两个临时变量是否能加快运行速度，这个例子体现不岀来，而下一个例子，则会有显著区别Plot(x,y)9由以上两个图形明显看出欧拉算法的精度远远不及龙格库塔法。前面的欧拉法只算了N=2F6个点，图像已经成了椭圆环了，而后面的龙格—库塔法，算了N=2F9，是前者的8倍，图像仍然是一个漂亮的椭圆"x=siny+cosxy=xcosy3、程序的优化F面我们用龙格-库塔法来画下面方程的轨线:不良编程的典范tich=0.001;N=2A19;t=O:h:N*h;x=zeros(size(t));y=zeros(size(t));x(1)=1；y(1)=0；fori=1:Nwx1=sin(y(i))+cos(x(i));wy1=x(i)*cos(y(i));wx2=sin(y(i)+wy1*h/2)+cos(x(i)+wx1*h/2);wy2=(x(i)+wx1*h/2)*cos(y(i)+wy1*h/2);wx3=sin(y(i)+wy2*h/2)+cos(x(i)+wx2*h/2);wy3=(x(i)+wx2*h/2)*cos(y(i)+wy2*h/2);wx4=sin(y(i)+wy3*h)+cos(x(i)+wx3*h);wy4=(x(i)+wx3*h)*cos(y(i)+wy3*h);x(i+1)=x(i)+(wx1+2*wx2+2*wx3+wx4)*h/6;y(i+1)=y(i)+(wy1+2*wy2+2*wy3+wy4)*h/6;end;Plot(x,y)toeElapsedtimeis15.092000seconds.1.6,,,,,TOC\o"1-5"\h\zIII1.4.-1.2.-1--0.8--0.6・・j0.4--0.2--f01L1111.522.533.5常量的使用tich=0.001;h2=h/2;h6=h/6;N=2A19;t=O:h:N*h;x=zeros(size(t));y=zeros(size(t));x(1)=1；y(1)=0；fori=1:Nwy仁x(i)*cos(y(i));wy2=(x(i)+wx1*h2)*cos(y(i)+wy1*h2);wy3=(x(i)+wx2*h2)*cos(y(i)+wy2*h2);wy4=(x(i)+wx3*h)*cos(y(i)+wy3*h);wx1=sin(y(i))+cos(x(i));wx2=sin(y(i)+wy1*h2)+cos(x(i)+wx1*h2);wx3=sin(y(i)+wy2*h2)+cos(x(i)+wx2*h2);wx4=sin(y(i)+wy3*h)+cos(x(i)+wx3*h);x(i+1)=x(i)+(wx1+2*wx2+2*wx3+wx4)*h6;y(i+1)=y(i)+(wy1+2*wy2+2*wy3+wy4)*h6;end;Plot(x,y)tocElapsedtimeis13.559000seconds.1.61.41.210.80.60.40.2011.522.533.5秒。现在我们看到，常量h2和h6的引入使得我们的一条轨道的计算时间减少了1.5我们还可以进一步改进，引进局部变量。使用局部变量tich=0.001;h2=h/2;h6=h/6;N=2A19;t=0:h:N*h;x=zeros(size(t));y=zeros(size(t));x(1)=1;y(1)=0;fori=1:Nux=x(i);ux=x(i)+wx1*h2;uy=y(i);wx1=sin(uy)+cos(ux);uy=y(i)+wy1*h2;wx2=sin(uy)+cos(ux);ux=x(i)+wx2*h2;uy=y(i)+wy2*h2;wx3=sin(uy)+cos(ux);ux=x(i)+wx3*h;uy=y(i)+wy3*h;wx4=sin(uy)+cos(ux);wy1=ux*cos(uy);wy2=ux*cos(uy);wy3=ux*cos(uy);wy4=ux*cos(uy);x(i+1)=x(i)+(wx1+2*wx2+2*wx3+wx4)*h6;y(i+1)=y(i)+(wy1+2*wy2+2*wy3+wy4)*h6;end;Plot(x,y)tocElapsedtimeis11.958000seconds.1.61.41.210.80.60.40.201.522.533.5(i)时间上的节省：我们看到，时间又节省了计算一条轨道，如果我们要画出一个参数为x'：次’x‘X3二Fcost'x=yy=-C(y_吸_妝3+FCOS^1t1.6秒，我们前后共节约时间3.1秒，不要小看这3.1秒，这只是[0,1；0,1]的舌头，如果取步长为0.001(这个不算很精确)，且不管计算Poincare截面所增加的大量时间，就算一个网格点算一条轨道，那么就要增加1000X1000X3.1/86400=35.8796(天！！！)可移植性：我们看到引入临时变量后，后面的这个方程和前面的例子改动的地方很少，而且极有规律这个例子中，临时变量的引入，之所以能使计算的时间缩短，在于龙格一库塔法中，现在的例子里每个ux和uy都要使用两次，因此还有一个原则：如果在循环体中某个与循环变量有关的量会用到不止一次，那就应该增加一个临时变量。道理弄明白了，那就只看大家的发挥了。四、一个具体的系统下面讲一个具体的系统：强迫Duffing方程，其中的程序由上面的讨论可以知道，可移植性很好，大家可以自由使用，但请不要改动其中的版权说明部分(一般发布源程序的人，都要写上这句话的！)，程序因为目前是在Word里运行，速度比较慢，所以运行长度都取得比较小，大家最好直接在MatLab里设置大一点的计算长度调用执行。这里的几个程序都是写成函数直接调用的，在MatLab里面的调用格式，和在这里用的是一样的。另外我不保证程序的正确性，那是你们的事情，切记切记！！！1、轨线图：Duffing(x0,y0,a,b,c,F,w,N,Dir)x0,y0显然为初始值，a,b,c,w分别对应于参数N为计算点数，Dir为方向。a=0.3;b=1;c=1;w=1.2;N=2A17;F=0.2;subplot(2,2,1);Duffing(-1.1,0,a,b,c,F,w,N);title('F=0.2');F=0.27;subplot(2,2,2);Duffing(-1.1,0,a,b,c,F,w,N);title('F=0.27');F=0.29;subplot(2,2,3);Duffing(-1.1,0,a,b,c,F,w,N);title('F=0.29');F=0.32;subplot(2,2,4);Duffing(-1.1,0,a,b,c,F,w,N);title('F=0.32');10.50-0.510.50-0.5F=0.2-0.5F=0.2710.50-0.5-1-1.5-1-0.5F=0.29F=0.32-0.52、变步长龙格库塔法轨线图：BBCDuffing(x0,y0,a,b,c,F,w,N)x0,y0显然为初始值，a,b,c,w分别对应于参数a=0.3;b=1;c=1;w=1.2;N=2A17;F=0.2;subplot(2,2,1);BBCDuffing(-1.1,0,a,b,c,F,w,N);title('F=0.2');F=0.27;subplot(2,2,2);BBCDuffing(-1.1,0,a,b,c,F,w,N);title('F=0.27');F=0.29;subplot(2,2,3);BBCDuffing(-1.1,0,a,b,c,F,w,N);title('F=0.29');■，N为次数。F=0.32;subplot(2,2,4);BBCDuffing(-1.1,0,a,b,c,F,w,N);title('F=0.32');-0.5F=0.2-0.5F=0.27-0.5-1—-1.5-1-0.510.50F=0.29-0.510.50-0.5-1-2-1012F=0.32变步长的效杲从前三个图中看不出来，因为，前三个图是周期的，而第四个图就有明显的差别了，我们作了同样的迭代次数，但是BBC的轨道显然走了更长的距离。3、Poincare映射：DuffingFlash(x0,y0,a,b,c,F,w,Delta,Count)其它参数同上，Delta为频闪步长，当Delta等于驱动外力周期时，就是Poincare截面，Count为频闪取点数，计算步长是程序自己内部自适应调整的，无需设置。a=0.3;b=1;c=1;w=1.2;Count=100;F=0.2;Delta=2*pi/w;subplot(2,2,1);DuffingFlash(-1.1,0,a,b,c,F,w,Delta,Count);F=0.27;Delta=2*pi/w;subplot(2,2,2);DuffingFlash(-1.1,0,a,b,c,F,w,Delta,Count);F=0.29;Delta=2*pi/w;subplot(2,2,3);DuffingFlash(-1.1,0,a,b,c,F,w,Delta,Count);F=0.32;Delta=2*pi/w;subplot(2,2,4);DuffingFlash(-1.1,0,a,b,c,F,w,Delta,500);0.630.70.6250.620.6150.61-0.9-0.8-0.7-0.60.605-10.20.40.60.810.60.50.4为确定前面三种情况确实是周期函数，上面的Poincare映射图给了我们清楚的图像。注意第一个图其实只有一个点。3、功率谱：DuffingGLP(x0,y0,a,b,c,F,w,N,Count)N是计算功率谱的长度，尽可能设置成2的幕，这是由快速傅立叶变换决定的。Count为要显示的功率谱的点数。程序是调用快速傅立叶变换，其实MatLab里有专门的计算功率谱的函数PSD还有些别的函数，因为最近没有时间去试，有时间的人自己去看MatLab的在线帮助。a=0.3;b=1;c=1;w=1.2;N=2A18;Count=100;F=0.2;subplot(2,2,1);DuffingGLP(-1.1,0,a,b,c,F,w,N,Count);F=0.27;subplot(2,2,2);DuffingGLP(-1.1,0,a,b,c,F,w,N,Count);F=0.29;subplot(2,2,3);DuffingGLP(-1.1,0,a,b,c,F,w,N,Count);F=0.32;subplot(2,2,4);DuffingGLP(-1.1,0,a,b,c,F,w,N,Count);功率谱图表明，解不断出现新的分频部分，最后呈现没有周期的运动。系统从什么时候开始轨线变得不是周期的呢？这就用到我们的分叉图。4、分叉图：DuffingFCH(xO,yO,a,b,c,w,FMin,FMax,dF,Count)FMin,FMaxdF分别为可调整参数的最小值，最大值和计算步长。这三个值决定了横向所要显示的点数，Count为纵向所要显示的点数，也就是Poincare截面的点数。目前下面的图是有问题的，FvO.2时，有好几条线，那是因为前面写程序的时候，暂态过程忽略过少，现在程序内部已经增加了暂态过程忽略的长度，但是这个程序运行太费时间，所以请大家自己去运行，我就不再运行一次了。a=0.3;b=1;c=1;w=1.2;FMin=0;FMax=0.4;dF=0.001;Count=100;DuffingFCH(-1.2,0,a,b,c,w,FMin,FMax,dF,Count);1鋼Tihiviidii3^ill从图中可以看到，系统大概在F=0.3的地方开始变得没有周期了。5、Liapunov指数DuffingLZS(x0,y0,a,b,c,w,F,Count);因为程序当时写的时候只为教学用的，所以除了系统的几个参数外，只设置了一个可调整的参数Count，Count是以tao作为步长，所要计算的点数，这里的tao是一个常数，知道李指数计算的人就知道这个数的意义，它本来也应该可以作为参数调整的，目前程序里tao=1另外本来还应该有个可调整参数，就是暂态过程所需要忽略的点数，目前程序里设置为200000,需要调整的人，可以将这两个参数也写成调用参数就行了。这个是很简单的工作，就留给没有写过函数的同仁动手试验吧。a=0.3;b=1;c=1;w=1.2;F=0.32;Count=4000;DuffingLZS(-1.1,0,a,b,c,w,F,Count);E=0.1306-0.43060此时Liapunov指数为(+,0,—)，因此确为混沌运动。6、轨线图：BBCDuffing(x0,y0,a,b,c,F,w,N)x0,y0显然为初始值，a,b,c,w分别对应于参数I,、.，」，，，N为ci次数。a=0.3;b=1;c=1;w=1.2;N=2A17;F=0.2;subplot(2,2,1);Duffing(-1.1,0,a,b,c,F,w,N);title('F=0.2');%F=0.27;subplot(2,2,2);Duffing(-1.1,0,a,b,c,F,w,N);title('F=0.27');%F=0.29;subplot(2,2,3);Duffing(-1.1,0,a,b,c,F,w,N);title('F=0.29');%F=0.32;subplot(2,2,4);Duffing(-1.1,0,a,b,c,F,w,N);title('F=0.32');