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

MMaattLLaabb编编程程简简介介一、先介绍MatLab里几个挺好用，但是大家以前可能没用到的命令，用过MatLab的从例子都应该能看明白用法，所以我基本不加注释。1、求解初等方程：022＋solve('x^2+w^2')ans=i*w-i*w2、求解微分方程：xx2,以及求函数的导数。xt=dsolve('D2x=-w^2*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、符号作图：)sin(tAysymst;%这个命令用来声明下面要用到的符号变量。ezplot(sin(t+pi/6),[0,2*pi]);axis([0,2*pi,-1,1]);目录第2页（共16页）0123456-1-0.8-0.6-0.4-0.200.20.40.60.81tsin(t+1/6)4、符号作带等高线曲面图：2222121xyzsymsxyw;ezsurfc(y*y/2+2*x*x,[-4,4,-8,8]);5、画向量场xyxQyyxPsin),(,),(u=-4:0.3:4;非线性动力系统第3页（共16页）[x,y]=meshgrid(u);%这个命令用来先生成一个数据网格以画向量场。quiver(x,y,y,-sin(x));-5-4-3-2-1012345-5-4-3-2-1012346、极坐标作图：211erTheta=0:0.01:5*pi;polar(Theta,1./sqrt(1+exp(-2*Theta)),'r');0.20.40.60.8130210602409027012030015033018007、画等高线：222)1(yxxz目录第4页（共16页）symsxy;ezcontour(x*x*(x-1)*(x-1)+y*y,[-0.3,1.3,-0.5,0.5]);xyx2(x-1)2+y2-0.200.20.40.60.811.2-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5二、MatLab的一些目前可能用到的命令列表，需要知道更详细的格式，可以在MatLab中查在线帮助。ezmesh画网格图ezmeshc画带等高线的网格图ezsurf画曲面图ezsurfc画带等高线的曲面图ezplot符号曲线图ezplot3三维符号曲线图ezpolar极坐标曲线图ezcontour画等高线图ezcontourf画填充等高线图plot数值二维图poar数值极坐标图plot3数值三维图mesh数值网格图meshc数值带等高线图surf数值曲面图surfc数值带等高线图quiver矢量图solve符号解代数方程dsolve符号解微分方程meshgrid产生数据网格gradient求数值梯度非线性动力系统第5页（共16页）三、微分方程的数值解法微分方程),(utfu，其中u为矢量。以下以简谐振动方程为例：xyyx4记)(mmtuu1、欧拉法：hutfuummmm),(1(1)h=0.001;N=2^16;t=0:h:N*h;x=zeros(size(t));y=zeros(size(t));%可以直接取x=zeros(1,N+1);但是一来容易范少算一个点的错误，二来下次改t的长%度时,还得记得改x的长度，这样编程x的长度随时间变量t自动取比较好x(1)=1;y(1)=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)目录第6页（共16页）-1.5-1-0.500.511.5-2.5-2-1.5-1-0.500.511.522.52、龙格库塔法：)22(6143211mmmmmmhuu(B.2)),()21,21()21,21(),(3423121huhtfhuhtfhuhtfutfmmmmmmmmmmmmmmm(B.3)h=0.001;h2=h/2;h6=h/6;N=2^19;%此处设置两个变量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;非线性动力系统第7页（共16页）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)-1-0.8-0.6-0.4-0.200.20.40.60.81-2-1.5-1-0.500.511.52由以上两个图形明显看出欧拉算法的精度远远不及龙格库塔法。前面的欧拉法只算了N=2^16个点，图像已经成了椭圆环了，而后面的龙格－库塔法，算了N=2^19，是前者的8倍，图像仍然是一个漂亮的椭圆3、程序的优化下面我们用龙格－库塔法来画下面方程的轨线：yxyxyxcoscossin(1)不良编程的典范tich=0.001;N=2^19;t=0: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;目录第8页（共16页）y(i+1)=y(i)+(wy1+2*wy2+2*wy3+wy4)*h/6;end;plot(x,y)tocElapsedtimeis15.092000seconds.11.522.533.500.20.40.60.811.21.41.6(2)常量的使用tich=0.001;h2=h/2;h6=h/6;N=2^19;t=0: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*h2)+cos(x(i)+wx1*h2);wy2=(x(i)+wx1*h2)*cos(y(i)+wy1*h2);wx3=sin(y(i)+wy2*h2)+cos(x(i)+wx2*h2);wy3=(x(i)+wx2*h2)*cos(y(i)+wy2*h2);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)*h6;y(i+1)=y(i)+(wy1+2*wy2+2*wy3+wy4)*h6;end;plot(x,y)tocElapsedtimeis13.559000seconds.非线性动力系统第9页（共16页）11.522.533.500.20.40.60.811.21.41.6现在我们看到，常量h2和h6的引入使得我们的一条轨道的计算时间减少了1.5秒。我们还可以进一步改进，引进局部变量。(3)使用局部变量tich=0.001;h2=h/2;h6=h/6;N=2^19;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=sin(uy)+cos(ux);wy1=ux*cos(uy);ux=x(i)+wx1*h2;uy=y(i)+wy1*h2;wx2=sin(uy)+cos(ux);wy2=ux*cos(uy);ux=x(i)+wx2*h2;uy=y(i)+wy2*h2;wx3=sin(uy)+cos(ux);wy3=ux*cos(uy);ux=x(i)+wx3*h;uy=y(i)+wy3*h;wx4=sin(uy)+cos(ux);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.目录第10页（共16页）11.522.533.500.20.40.60.811.21.41.6(i)时间上的节省：我们看到，时间又节省了1.6秒，我们前后共节约时间3.1秒，不要小看这3.1秒，这只是计算一条轨道，如果我们要画出一个参数为]1,0;1,0[的舌头，如果取步长为0.001（这个不算很精确），且不管计算Poincare截面所增加的大量时间，就算一个网格点算一条轨道，那么就要增加1000×1000×3.1/86400＝35.8796(天！！！)(ii)可移植性：我们看到引入临时变量后，后面的这个方程和前面的例子改动的地方很少，而且极有规律(iii)这个例子中，临时变量的引入，之所以能使计算的时间缩短，在于龙格－库塔法中，现在的例子里每个ux和uy都要使用两次，因此还有一个原则：如果在循环体中某个与循环变量有关的量会用到不止一次，那就应该增加一个临时变量。道理弄明白了，那就只看大家的发挥了。四、一个具体的系统下面讲一个具体的系统：强迫Duffing方程，其中的程序由上面的讨论可以知道，可移植性很好，大家可以自由使用，但请不要改动其中的版权说明部分(一般发布源程序的人，都要写上这句话的！)，程序因为目前是在Word里运行，速度比较慢，所以运行长度都取得比较小，大家最好直接在MatLab里设置大一点的计算长度调用执行。这里的几个程序都是写成函数直接调用的，在MatLab里面的调用格式，和在这里用的是一样的。另外我不保证程序的正确性，那是你们的事情，切记切记！！！tFxxxxcos3或tFxxyyyxcos31、轨线图：Duffing(x0,y0,a,b,c,F,w,N,Dir)x0,y0显然为初始值，a,b,c,w分别对应于参数，，，，N为计算点数，Dir为方向。非线性动力系统第11页（共16页）a=0.3;b=1;c=1;w=1.2;N=2^17;F=0.2;subplot(2,2,1);Duffing(-1.1,0,a,b,c