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

1、精选优质文档-倾情为你奉上MatLab编程简介一、 先介绍MatLab里几个挺好用，但是大家以前可能没用到的命令，用过MatLab的从例子都应该能看明白用法，所以我基本不加注释。1、 求解初等方程：solve(x2+w2) ans = i*w -i*w 2、 求解微分方程：,以及求函数的导数。xt=dsolve(D2x = -w2*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)*w 3、符号作图：syms t;%这个命令用来声明下面要用到的符号变量。ezplot(sin(t+pi/6),0,

2、2*pi);axis(0,2*pi,-1,1); 4、符号作带等高线曲面图：syms x y w; ezsurfc(y*y/2+2*x*x,-4,4,-8,8); 5、画向量场u=-4:0.3:4;x,y=meshgrid(u);%这个命令用来先生成一个数据网格以画向量场。quiver(x,y,y,-sin(x); 6、极坐标作图：Theta=0:0.01:5*pi; polar(Theta,1./sqrt(1+exp(-2*Theta),r); 7、画等高线： syms x y;ezcontour(x*x*(x-1)*(x-1)+y*y,-0.3,1.3,-0.5,0.5); 二、MatLa

3、b的一些目前可能用到的命令列表，需要知道更详细的格式，可以在MatLab中查在线帮助。ezmesh画网格图ezmeshc画带等高线的网格图ezsurf画曲面图ezsurfc画带等高线的曲面图ezplot符号曲线图ezplot3三维符号曲线图ezpolar极坐标曲线图ezcontour画等高线图ezcontourf画填充等高线图plot数值二维图poar数值极坐标图plot3数值三维图mesh数值网格图meshc数值带等高线图surf数值曲面图surfc数值带等高线图quiver矢量图solve符号解代数方程dsolve符号解微分方程meshgrid产生数据网格gradient求数值梯度三、微分

4、方程的数值解法微分方程 ，其中为矢量。以下以简谐振动方程为例：记1、欧拉法： (1)h=0.001;N=216;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;for i=1:N %这里的循环次数极容易弄错，要么多算要么少算一个点，这种错误往往还查不出来。x(i+1)=x(i)+y(i)*h;y(i+1)=y(i)-4*x(i)*h;end;plot(x,y) 2、龙格库塔法

5、： (B.2) (B.3)h=0.001;h2=h/2;h6=h/6;N=219;%此处设置两个变量h2和h6是为了减少下面的大循环中的除法的次数，记住一个原则：如果在循环里有某个量%与循环变量无关,那么所有关于这个量的运算能算好的都先算好。t=0:h:N*h;x=zeros(size(t);y=zeros(size(t);x(1)=1;y(1)=0;for i=1:N ux=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; ux=x(i)+wx2*h2; uy=y(i)+wy

6、2*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,

7、y) 由以上两个图形明显看出欧拉算法 的精度远远不及龙格库塔法。前面的欧拉法只算了N=216个点，图像已经成了椭圆环了，而后面的龙格库塔法，算了N=219，是前者的8倍，图像仍然是一个漂亮的椭圆3、程序的优化下面我们用龙格库塔法来画下面方程的轨线：(1) 不良编程的典范tich=0.001;N=219;t=0:h:N*h;x=zeros(size(t);y=zeros(size(t);x(1)=1;y(1)=0;for i=1:N wx1=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);

8、 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) toc Elapsed time is

9、15. seconds. (2) 常量的使用tich=0.001;h2=h/2;h6=h/6;N=219;t=0:h:N*h;x=zeros(size(t);y=zeros(size(t);x(1)=1;y(1)=0;for i=1:N wx1=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)+w

10、y2*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) toc Elapsed time is 13. seconds. 现在我们看到，常量h2和h6的引入使得我们的一条轨道的计算时间减少了1.5秒。我们还可以进一步改进，引进局部变量。(3)使用局部变量tich=0.001;h2=h/2;h6=h/6;N=219;t=0

11、:h:N*h;x=zeros(size(t);y=zeros(size(t);x(1)=1;y(1)=0;for i=1:N ux=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

12、);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) toc Elapsed time is 11. seconds. (i)时间上的节省：我们看到，时间又节省了1.6秒，我们前后共节约时间3.1秒，不要小看这3.1秒，这只是计算一条轨道，如果我们要画出一个参数为的舌头，如果取步长为0.001（这个不算很精确），且不管计算Poincare截面所增加的大量时间，就算一个网格点算一条轨道，那么就要增加100010003.1/8640035.

13、8796(天！)(ii)可移植性：我们看到引入临时变量后，后面的这个方程和前面的例子改动的地方很少，而且极有规律(iii) 这个例子中，临时变量的引入，之所以能使计算的时间缩短，在于龙格库塔法中，现在的例子里每个ux和uy都要使用两次，因此还有一个原则：如果在循环体中某个与循环变量有关的量会用到不止一次，那就应该增加一个临时变量。道理弄明白了，那就只看大家的发挥了。四、一个具体的系统下面讲一个具体的系统：强迫Duffing方程，其中的程序由上面的讨论可以知道，可移植性很好，大家可以自由使用，但请不要改动其中的版权说明部分(一般发布源程序的人，都要写上这句话的！)，程序因为目前是在Word里运行

14、，速度比较慢，所以运行长度都取得比较小，大家最好直接在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=217;F=0.2; subplot(2,2,1);Duffing(-1.1,0,a,b,c,F,w,N);title(F=0.2);F=0.27;su

15、bplot(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); 2、 变步长龙格库塔法轨线图：BBCDuffing(x0,y0,a,b,c,F,w,N) x0,y0显然为初始值，a,b,c,w分别对应于参数，N为次数。a=0.3;b=1;c=1;w=1.2;N=217;F=0.2; subplo

16、t(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);F=0.32;subplot(2,2,4);BBCDuffing(-1.1,0,a,b,c,F,w,N);title(F=0.32); 变步长的效果从前三个图中看不出来，因为，前三个图是周期的，而第四个图就有明显的差别了，我们

17、作了同样的迭代次数，但是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;subplo

18、t(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); 为确定前面三种情况确实是周期函数，上面的Poincare映射图给了我们清楚的图像。注意第一个图其实只有一个点。3、 功率谱：DuffingGLP(x0,y0,a,b,c,F,w,N,

19、Count)N是计算功率谱的长度，尽可能设置成2的幂，这是由快速傅立叶变换决定的。Count为要显示的功率谱的点数。程序是调用快速傅立叶变换，其实MatLab里有专门的计算功率谱的函数PSD，还有些别的函数，因为最近没有时间去试，有时间的人自己去看MatLab的在线帮助。a=0.3;b=1;c=1;w=1.2;N=218;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;sub

20、plot(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(x0,y0,a,b,c,w,FMin,FMax,dF,Count)FMin，FMax，dF分别为可调整参数的最小值，最大值和计算步长。这三个值决定了横向所要显示的点数，Count为纵向所要显示的点数，也就是Poincare截面的点数。目前下面的图是有问题的，F0.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); 从图中可以看到，系统大概在的地方开始变得没有周期了。5、 Liapunov指数DuffingLZS(x0,y0,a,b,c,w,F,Count);