非线性动力学作业.doc

非线性动力学第一次上机大作业姓名：王建春学号：p2014020471、 题目分析在平衡点的稳定性其中方程是1，u=0.1,w0=0.22, u2=w02=0.013, u=0.02,w0=0.01初值是x1=(0.1,0.5,1.5),x2=(0.2,0.6,1.8)2、 分析 转换格式为如下形式 3、 利用MATLAB编程求解程序如下1，创建函数xprim2，并将其保存在M文件xprim2.m中function xprim = xprim2(t,x)xprim = x(2);-2*u*x(2)-w2*x(1);2，然后调用一个求特征值与ODE算法和画出解的图形symsx1 x2;f=x2;-2*u*x2-w2*x1;v=x1,x2;jacob=jacobian(f,v);a=eig(jacob);t,x=ode45(xprim2,0 200,a;b)Plot(t,x);xlabel(time t0=0,tt=200);ylabel(x value x1(0)=a,x2(0)=b);Plot(x(:,2),x(:,1);注：对于不同的初始条件和不同的控制方程需要带入不同的u,w,a,b即可求得想要的结果。4、 试验结果1、 在M文件中以u=0.1,w=0.2得到三个不同初值x1=(0.1,0.5,1.5)和x2=(0.2,0.6,1.8)的结果如下图，其特征值是a=-1/10+1/10*i*3(1/2)和-1/10-1/10*i*3(1/2)初值x1=0.1和x2=0.2随时间变化的曲线图初值x1=0.1和x2=0.2的相位曲线图初值x1=0.5，x2=0.6随时间变化的曲线图初值x1=0.5,x2=0.6的相位曲线图初值x1=1.5,x2=1.8随时间变化的曲线图初值x1=1.5,x2=1.8的相位曲线图2、 在M文件中以u=0.02,w=0.01,得到三个不同初值x1=(0.1,0.5,1.5)和x2=(0.2,0.6,1.8)的结果如下图,其特征值是a=-1/50+1/100*3(1/2)和-1/50-1/100*3(1/2)。初值x1=0.1和x2=0.2随时间变化的曲线图初值x1=0.1和x2=0.2的相位曲线图初值x1=0.5和x2=0.6随时间变化的曲线图初值x1=0.5和x2=0.6的相位曲线图初值x1=1.5和x2=1.8随时间变化的曲线图初值x1=1.5和x2=1.8的相位曲线图3、 在M文件中以u=0.1,w2=0.01得到三个不同初值x1=(0.1,0.5,1.5)和x2=(0.2,0.6,1.8)的结果如下图，其特征值是-1/10不同初值所得到的x1,x2随时间变化曲线不同初值的相位曲线图4、 在M文件中以u=-0.1,w2=0.01得到三个不同初值x1=(0.1,0.5,1.5)和x2=(0.2,0.6,1.8)的结果如下图特征值是1/10.不同初值的x1，x2随时间变化的曲线图x1，x2的相位曲线图五、结果分析首先令方程组右边=0得到平衡点是（0,0）第一种情况下，u=0.1,w=0.2求得jacobian的特征值是-1/10+1/10*i*3(1/2)和-1/10-1/10*i*3(1/2)由于实部是负值，w0,此时平衡点为稳定的焦点，由上面图形曲线可以得出，平衡点也可以称为汇。第二种情况下，u=0.02,w=0.01，求得jacobian的特征值是-1/50+1/100*3(1/2)和-1/50-1/100*3(1/2)，都是实数平衡点是稳定的双曲结点。第三种情况下，u2=w2=0.01,取u=0.1和u=-0.1下分别得到特征值是-1/10和1/10，得出平衡点分别是稳定的结点和不稳点的结点。#include #include #define N 3#define h 1/(N+1) /*步长*/#define M 50main() double jacobi(int AN*NN*N,double bN*N);int i,j,L,l,r,up,d,bj=0;double uN+2N+2=0; /*包含边界点*/double vN*N;int AN*NN*N=0;double BN*N;/*初始化矩阵B*/for(i=0;iN*N;i+) Bi=0.25*0.25; /*按列编号*/ for(i=1;i=N;i+) for(j=1;j=N;j+) for(L=0;LN*N;L+) vL=uij; /*给矩阵A赋值*/ for(L=0;LN*N;L+) ALL=4; if(L/N=0) BL=BL+bj; else l=N*(L/N-1)+L%N; ALl=-1; if(L/N+2=N) BL=BL+bj; else r=N*(L/N+1)+L%N; ALr=-1; if(L%N+2=N) BL=BL+bj; elseup=N*(L/N)+L%N+1; ALup=-1; if(L%N=0) BL=BL+bj; else d=N*(L/N)+L%N-1; ALd=-1; /*输出矩阵A*/ printf(系数矩阵A: );printf(n); for(i=0;iN*N;i+) for(j=0;jN*N;j+) printf(%5d ,Aij); printf(n); printf(n);/*输出矩阵b*/printf(输出b: );printf(n);for(i=0;iN*N;i+)printf(%f ,Bi);printf(n);/*调用函数计算结果*/ vN*N=jacobi(A,B); system(pause);double jacobi(int AN*NN*N,double bN*N) double DN*NN*N=0;double LN*NN*N=0;double UN*NN*N=0;double BN*NN*N=0;int i,j,m;double gN*N=0;double xMN*N=1;double vN*N;/*给D矩阵赋值*/for(i=0;iN*N;i+) /*Dii=1/Aii;*/ Dii=0.25;/*输出矩阵Dprintf(输出矩阵D: );printf(n); for(i=0;iN*N;i+) for(j=0;jN*N;j+) printf(%f ,Dij); printf(n); printf(n); */*给L矩阵赋值*/for(i=0;iN*N;i+) for(j=0;ji;j+) Lij=-Aij; /*输出L矩阵printf(输出矩阵L: );printf(n); for(i=0;iN*N;i+) for(j=0;jN*N;j+) printf(%f ,Lij); printf(n); printf(n); */*给U矩阵赋值*/for(j=0;jN*N;j+) for(i=0;ij;i+) Uij=-Aij; /*输出U矩阵printf(输出矩阵U: );printf(n); for(i=0;iN*N;i+) for(j=0;jN*N;j+) printf(%f ,Uij); printf(n); printf(n); */ /*给B矩阵赋值*/for(i=0;iN*N;i+) for(m=0;mN*N;m+) for(j=0;jN*N;j+) Bim+=Dij*(Ljm+Ujm); /*输出B矩阵printf(输出矩阵B: );printf(n); for(i=0;iN*N;i+) for(j=0;jN*N;j+) printf(%f ,Bij); printf(n); printf(n); */*给g矩阵赋值*/for(i=0;iN*N;i+) for(j=0;jN*N;j+) gi+=Dij*bj; /*输出g矩阵printf(输出矩阵g: );printf(n); for(i=0;iN*N;i+) printf(%f,gi); printf(n); printf(n); */*jacobi迭代格式*/for(i=0;iM;i+) for(m=0;mN*N;m+) for(j=0;jN*N;j+) xi+1m+=Bmj*xij; xi+1m+=gm; /*输出X矩阵printf(输出矩阵X: );printf(n); for(i=0;iM;i+) for(j=0;j