谁有计算4维混沌系统lyapunov指数的程序?
如<BR>dx/dt=-y-w;<BR>dy/dt=x+a*y+z;<BR>dz/dt=b*z-c*w;<BR>dw/dt=x*w+d;<BR><BR>a = 0.25;<BR>b = 0.05;<BR>c = 0.5;<BR>d = 3;不知道
<STRONG>你有计算３维混沌系统lyapunov指数的程序吗?</STRONG>largest lyapunov exponent duffing system
<P>#include <iostream.h><BR>#include <fstream.h><BR>#include <math.h><BR>#defineM 1000<BR>#defineN1 50<BR>#defineN2 50<BR>#definePI 3.14159265</P><P>double *f(double t,double x)<BR>{<BR> double a=0.13,e=1.0,w=1.0,y=0.18,d=0.02;<BR> double *p,m;<BR> m=x;<BR> m=x-a*x-e*x*(x*x+x*x)+y*cos(w*t);<BR> m=x;<BR> m=x-a*x-e*x*(x*x+x*x)+y*cos(w*t);<BR> p=m;<BR> return p;<BR>}</P>
<P>double *Runge_Kutta(double t0,double t1,double x00,double x01,double x02,double x03,double h)<BR>{<BR> //Runge_Kutta法解微分方程组<BR> //x'=f(x) length(x)=4<BR>double t,x,*p;<BR>x=x00;<BR>x=x01;<BR>x=x02;<BR>x=x03;<BR>double *k1,*k2,*k3,*k4;<BR>for(t=t0;t<t1;t=t+h)<BR>{<BR> k1=f(t,x);<BR> x+=h*k1/2;<BR> x+=h*k1/2;<BR>x+=h*k1/2;<BR>x+=h*k1/2;<BR> k2=f(t+h/2,x);<BR> x+=h*k2/2;<BR> x+=h*k2/2;<BR>x+=h*k2/2;<BR>x+=h*k2/2;<BR> k3=f(t+h/2,x);<BR> x+=h*k3;<BR> x+=h*k3;<BR>x+=h*k3;<BR>x+=h*k3;<BR> k4=f(t+h,x);<BR> x+=h*(k1+2*k2+2*k3+k4)/6;<BR> x+=h*(k1+2*k2+2*k3+k4)/6;<BR>x+=h*(k1+2*k2+2*k3+k4)/6;<BR>x+=h*(k1+2*k2+2*k3+k4)/6;<BR>}<BR>p=x;<BR> return p;<BR>}<BR>double Lyapunov_Exponent(double x,double z)<BR>{<BR> double delt=2*PI;<BR> double d0=sqrt((x-z)*(x-z)+(x-z)*(x-z));<BR> double y,d;<BR> for(int i=0;i<M;i++)<BR> d=0;<BR> for(i=0;i<M;i++)<BR> {<BR> x=Runge_Kutta(i*delt,(i+1)*delt,x,x,x,x,0.001);<BR> x=Runge_Kutta(i*delt,(i+1)*delt,x,x,x,x,0.001);<BR> y=Runge_Kutta(i*delt,(i+1)*delt,z,z,z,z,0.001);<BR> y=Runge_Kutta(i*delt,(i+1)*delt,z,z,z,z,0.001);<BR> d=sqrt((x-y)*(x-y)+(x-y)*(x-y));<BR> z=x+(d0/d)*(x-y);<BR> z=x+(d0/d)*(x-y);<BR> }<BR> double le=0;<BR> for(i=0;i<M;i++)<BR> le+=log(d/d0);<BR> le/=(M*delt);<BR> return le;<BR>}<BR>void main()<BR>{<BR> float n=2,s=-2,w=-2,e=2;<BR> float h1=(e-w)/N1,h2=(n-s)/N2;<BR>double le;<BR>double x,z;<BR>for(int i=0;i<N1;i++)<BR>for(int j=0;j<N2;j++)<BR>le=0;</P>
<P>for(i=0;i<N1;i++)<BR>for(int j=0;j<N2;j++)<BR>{<BR> x=i*h1+w;<BR> x=i*h2+s;<BR> z=(i+1)*h1+w;<BR> z=(i+1)*h2+s;<BR> le=Lyapunov_Exponent(x,z);<BR> cout<<"i="<<i<<"\tj="<<j<<"\n";<BR>}<BR> ofstream file;<BR>file.open("Lyapunov_Exponent.txt");<BR>if(!file)<BR>{<BR> cout<<"can't open file Lyapunov_Exponent.txt\n";<BR>}<BR>for(i=0;i<N1;i++)<BR>{<BR>for(int j=0;j<N2;j++)<BR>{<BR> cout<<le<<"\t";<BR> file<<le;<BR> if(j==N2-1)<BR> file<<";";<BR> else<BR> file<<",";<BR>}<BR>cout<<"\n";<BR>}</P>
<P>}<BR></P> 有用Matlab语言编写的? matlab太慢了。。。。。。。。。。。。。 好像四维和二维的差别不大吧 把计算方法搞搞清楚，还是自己写一个小程序解决问题来的实在^_^ Duffing方程算3维的还是4维的 原帖由 ddb_2005 于 2006-8-21 23:43 发表
Duffing方程算3维的还是4维的
Duffing 振子、van der Pol 振子等非自治系统可以是二维系统 这段C语言的程序貌似有些错误
我修改了一下，不知道意思对不对，不过倒是可以运行了；
目前还没算出结果，算出来了再把结果贴出来……
就是程序没有说明，郁闷……
visual c++ 6.0
--------------------------------
// C_Lyapunov.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
// largest lyapunov exponent duffing system
#include <iostream.h>
#include <fstream.h>
#include <math.h>
#define M 1000
#define N1 50
#define N2 50
#define PI 3.14159265
double *f(double t,double x)
{
double a=0.13,e=1.0,w=1.0,y=0.18,d=0.02;
double *p,m;
m=x;
m=x-a*x-e*x*(x*x+x*x)+y*cos(w*t);
m=x;
m=x-a*x-e*x*(x*x+x*x)+y*cos(w*t);
p=m;
return p;
}
double *Runge_Kutta(double t0,double t1,double x00,double x01,double x02,double x03,double h)
{
//Runge_Kutta法解微分方程组
//x'=f(x) length(x)=4
double t,x,*p;
x=x00;
x=x01;
x=x02;
x=x03;
double *k1,*k2,*k3,*k4;
for(t=t0;t<t1;t=t+h)
{
k1=f(t,x);
x+=h*k1/2;
x+=h*k1/2;
x+=h*k1/2;
x+=h*k1/2;
k2=f(t+h/2,x);
x+=h*k2/2;
x+=h*k2/2;
x+=h*k2/2;
x+=h*k2/2;
k3=f(t+h/2,x);
x+=h*k3;
x+=h*k3;
x+=h*k3;
x+=h*k3;
k4=f(t+h,x);
x+=h*(k1+2*k2+2*k3+k4)/6;
x+=h*(k1+2*k2+2*k3+k4)/6;
x+=h*(k1+2*k2+2*k3+k4)/6;
x+=h*(k1+2*k2+2*k3+k4)/6;
}
p=x;
return p;
}
double Lyapunov_Exponent(double x,double z)
{
double delt=2*PI;
double d0=sqrt((x-z)*(x-z)+(x-z)*(x-z));
double y,d;
for(int i=0;i<M;i++)
d=0;
// double y, d;
//
// for(int i=0;i<M;i++)
// d = 0;
for(i=0;i<M;i++)
{
//cout<<"No."<<i<<endl;
x=Runge_Kutta(i*delt,(i+1)*delt,x,x,x,x,0.001);
x=Runge_Kutta(i*delt,(i+1)*delt,x,x,x,x,0.001);
y=Runge_Kutta(i*delt,(i+1)*delt,z,z,z,z,0.001);
y=Runge_Kutta(i*delt,(i+1)*delt,z,z,z,z,0.001);
d=sqrt((x-y)*(x-y)+(x-y)*(x-y));
z=x+(d0/d)*(x-y);
z=x+(d0/d)*(x-y);
}
double le=0;
for(i=0;i<M;i++)
le+=log(d/d0);
le/=(M*delt);
return le;
}
void main()
{
float n=2,s=-2,w=-2,e=2;
float h1=(e-w)/N1,h2=(n-s)/N2;
double le;
double x,z;
for(int i=0;i<N1;i++)
{
for(int j=0;j<N2;j++)
le=0;
}
for(i=0;i<N1;i++)
{
for(int j=0;j<N2;j++)
{
cout<<"row:"<<i<<'\t'<<"col:"<<j<<" : Lyapunov_Exponent"<<endl;
x=i*h1+w;
x=i*h2+s;
z=(i+1)*h1+w;
z=(i+1)*h2+s;
// cout<<"Lyapunov_Exponent"<<endl;
le=Lyapunov_Exponent(x,z);
cout<<"i="<<i<<"\tj="<<j<<"\n";
}
}
ofstream file;
file.open("Lyapunov_Exponent.txt");
if(!file)
{
cout<<"can't open file Lyapunov_Exponent.txt\n";
}
for(i=0;i<N1;i++)
{
for(int j=0;j<N2;j++)
{
cout<<le<<"\t";
file<<le;
if(j==N2-1)
file<<";";
else
file<<",";
}
cout<<"\n";
}
printf("Hello World!\n");
// return 0;
}
// int main(int argc, char* argv[])
// {
// printf("Hello World!\n");
// return 0;
// }
回复 10楼 epistemer 的帖子
计算了好久，终于出结果了，贴出来大家看看对不对数据：
绘图： 挺好的，算是基本算法了
回复 12楼 xinwilliam 的帖子
怎么都是负的？ 图是怎么画的？matlab吗？
页:
[1]