声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 4234|回复: 13

[分形与混沌] 谁有计算4维混沌系统lyapunov指数的程序?

[复制链接]
发表于 2006-6-17 18:13 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
如<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;

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2006-6-20 08:06 | 显示全部楼层

不知道

<STRONG>你有计算3维混沌系统lyapunov指数的程序吗?</STRONG>
发表于 2006-6-29 09:49 | 显示全部楼层

largest lyapunov exponent duffing system

<P>#include &lt;iostream.h&gt;<BR>#include &lt;fstream.h&gt;<BR>#include &lt;math.h&gt;<BR>#define  M    1000<BR>#define  N1   50<BR>#define  N2   50<BR>#define  PI   3.14159265</P>
<P>double *f(double t,double x[4])<BR>{<BR>   double a=0.13,e=1.0,w=1.0,y=0.18,d=0.02;<BR>   double *p,m[4];<BR>   m[0]=x[1];<BR>   m[1]=x[0]-a*x[1]-e*x[0]*(x[0]*x[0]+x[2]*x[2])+y*cos(w*t);<BR>   m[2]=x[3];<BR>   m[3]=x[2]-a*x[3]-e*x[2]*(x[0]*x[0]+x[2]*x[2])+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[4],*p;<BR>x[0]=x00;<BR>x[1]=x01;<BR>x[2]=x02;<BR>x[3]=x03;<BR>double *k1,*k2,*k3,*k4;<BR>  for(t=t0;t&lt;t1;t=t+h)<BR>  {<BR>    k1=f(t,x);<BR>        x[0]+=h*k1[0]/2;<BR>     x[1]+=h*k1[1]/2;<BR>  x[2]+=h*k1[2]/2;<BR>  x[3]+=h*k1[3]/2;<BR>    k2=f(t+h/2,x);<BR>        x[0]+=h*k2[0]/2;<BR>     x[1]+=h*k2[1]/2;<BR>  x[2]+=h*k2[2]/2;<BR>  x[3]+=h*k2[3]/2;<BR>    k3=f(t+h/2,x);<BR>     x[0]+=h*k3[0];<BR>     x[1]+=h*k3[1];<BR>  x[2]+=h*k3[2];<BR>  x[3]+=h*k3[3];<BR>    k4=f(t+h,x);<BR>     x[0]+=h*(k1[0]+2*k2[0]+2*k3[0]+k4[0])/6;<BR>     x[1]+=h*(k1[1]+2*k2[1]+2*k3[1]+k4[1])/6;<BR>  x[2]+=h*(k1[2]+2*k2[2]+2*k3[2]+k4[2])/6;<BR>  x[3]+=h*(k1[3]+2*k2[3]+2*k3[3]+k4[3])/6;  <BR>  }<BR>  p=x;<BR> return p;<BR>}<BR>double Lyapunov_Exponent(double x[2],double z[2])<BR>{<BR> double delt=2*PI;<BR>    double d0=sqrt((x[0]-z[0])*(x[0]-z[0])+(x[1]-z[1])*(x[1]-z[1]));<BR> double y[2],d[M];<BR> for(int i=0;i&lt;M;i++)<BR>        d=0;<BR> for(i=0;i&lt;M;i++)<BR> {<BR>      x[0]=Runge_Kutta(i*delt,(i+1)*delt,x[0],x[1],x[0],x[1],0.001)[0];<BR>      x[1]=Runge_Kutta(i*delt,(i+1)*delt,x[0],x[1],x[0],x[1],0.001)[1];<BR>      y[0]=Runge_Kutta(i*delt,(i+1)*delt,z[0],z[1],z[0],z[1],0.001)[0];<BR>      y[1]=Runge_Kutta(i*delt,(i+1)*delt,z[0],z[1],z[0],z[1],0.001)[1];<BR>   d=sqrt((x[0]-y[0])*(x[0]-y[0])+(x[1]-y[1])*(x[1]-y[1]));<BR>          z[0]=x[0]+(d0/d)*(x[0]-y[0]);<BR>          z[1]=x[1]+(d0/d)*(x[1]-y[1]);<BR> }<BR> double le=0;<BR> for(i=0;i&lt;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[N1][N2];<BR>  double x[2],z[2];<BR>  for(int i=0;i&lt;N1;i++)<BR>  for(int j=0;j&lt;N2;j++)<BR>  le[j]=0;</P>
<P>  for(i=0;i&lt;N1;i++)<BR>  for(int j=0;j&lt;N2;j++)<BR>  {<BR>       x[0]=i*h1+w;<BR>       x[1]=i*h2+s;<BR>       z[0]=(i+1)*h1+w;<BR>       z[1]=(i+1)*h2+s;  <BR>    le[j]=Lyapunov_Exponent(x,z);<BR>    cout&lt;&lt;"i="&lt;&lt;i&lt;&lt;"\tj="&lt;&lt;j&lt;&lt;"\n";<BR>  }<BR> ofstream file;<BR>  file.open("Lyapunov_Exponent.txt");<BR>  if(!file)<BR>  {<BR>   cout&lt;&lt;"can't open file Lyapunov_Exponent.txt\n";<BR>  }<BR>  for(i=0;i&lt;N1;i++)<BR>  {<BR>  for(int j=0;j&lt;N2;j++)<BR>  {<BR>   cout&lt;&lt;le[j]&lt;&lt;"\t";<BR>         file&lt;&lt;le[j];<BR>   if(j==N2-1)<BR>    file&lt;&lt;";";<BR>   else<BR>    file&lt;&lt;",";<BR>  }<BR>  cout&lt;&lt;"\n";<BR>  }</P>
<P>}<BR></P>
发表于 2006-6-29 20:23 | 显示全部楼层
有用Matlab语言编写的?
发表于 2006-6-30 09:40 | 显示全部楼层
matlab太慢了。。。。。。。。。。。。。
发表于 2006-7-13 16:46 | 显示全部楼层
好像四维和二维的差别不大吧
发表于 2006-7-13 22:04 | 显示全部楼层
把计算方法搞搞清楚,还是自己写一个小程序解决问题来的实在^_^
发表于 2006-8-21 23:43 | 显示全部楼层
Duffing方程算3维的还是4维的
发表于 2006-8-22 07:50 | 显示全部楼层
原帖由 ddb_2005 于 2006-8-21 23:43 发表
Duffing方程算3维的还是4维的


Duffing 振子、van der Pol 振子等非自治系统可以是二维系统
发表于 2010-1-28 23:36 | 显示全部楼层
这段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[4])
{
        double a=0.13,e=1.0,w=1.0,y=0.18,d=0.02;
        double *p,m[4];
        m[0]=x[1];
        m[1]=x[0]-a*x[1]-e*x[0]*(x[0]*x[0]+x[2]*x[2])+y*cos(w*t);
        m[2]=x[3];
        m[3]=x[2]-a*x[3]-e*x[2]*(x[0]*x[0]+x[2]*x[2])+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[4],*p;
        x[0]=x00;
        x[1]=x01;
        x[2]=x02;
        x[3]=x03;
        double *k1,*k2,*k3,*k4;
        for(t=t0;t<t1;t=t+h)
        {
                k1=f(t,x);
                x[0]+=h*k1[0]/2;
                x[1]+=h*k1[1]/2;
                x[2]+=h*k1[2]/2;
                x[3]+=h*k1[3]/2;
                k2=f(t+h/2,x);
                x[0]+=h*k2[0]/2;
                x[1]+=h*k2[1]/2;
                x[2]+=h*k2[2]/2;
                x[3]+=h*k2[3]/2;
                k3=f(t+h/2,x);
                x[0]+=h*k3[0];
                x[1]+=h*k3[1];
                x[2]+=h*k3[2];
                x[3]+=h*k3[3];
                k4=f(t+h,x);
                x[0]+=h*(k1[0]+2*k2[0]+2*k3[0]+k4[0])/6;
                x[1]+=h*(k1[1]+2*k2[1]+2*k3[1]+k4[1])/6;
                x[2]+=h*(k1[2]+2*k2[2]+2*k3[2]+k4[2])/6;
                x[3]+=h*(k1[3]+2*k2[3]+2*k3[3]+k4[3])/6;
        }
        p=x;
        return p;
}
double Lyapunov_Exponent(double x[2],double z[2])
{
        double delt=2*PI;
        double d0=sqrt((x[0]-z[0])*(x[0]-z[0])+(x[1]-z[1])*(x[1]-z[1]));
        double y[2],d[M];
       
        for(int i=0;i<M;i++)
                d=0;
       
        //         double y[2], d;
        //
        //         for(int i=0;i<M;i++)
        //                 d = 0;
       
        for(i=0;i<M;i++)
        {
                //cout<<"No."<<i<<endl;
                x[0]=Runge_Kutta(i*delt,(i+1)*delt,x[0],x[1],x[0],x[1],0.001)[0];
                x[1]=Runge_Kutta(i*delt,(i+1)*delt,x[0],x[1],x[0],x[1],0.001)[1];
                y[0]=Runge_Kutta(i*delt,(i+1)*delt,z[0],z[1],z[0],z[1],0.001)[0];
                y[1]=Runge_Kutta(i*delt,(i+1)*delt,z[0],z[1],z[0],z[1],0.001)[1];
                d=sqrt((x[0]-y[0])*(x[0]-y[0])+(x[1]-y[1])*(x[1]-y[1]));
                z[0]=x[0]+(d0/d)*(x[0]-y[0]);
                z[1]=x[1]+(d0/d)*(x[1]-y[1]);
        }
        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[N1][N2];
        double x[2],z[2];
        for(int i=0;i<N1;i++)
        {       
                for(int j=0;j<N2;j++)
                        le[j]=0;
        }
        for(i=0;i<N1;i++)
        {       
               
                for(int j=0;j<N2;j++)
                {
                        cout<<"row:"<<i<<'\t'<<"col:"<<j<<" : Lyapunov_Exponent"<<endl;
                        x[0]=i*h1+w;
                        x[1]=i*h2+s;
                        z[0]=(i+1)*h1+w;
                        z[1]=(i+1)*h2+s;
                //         cout<<"Lyapunov_Exponent"<<endl;
                        le[j]=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[j]<<"\t";
                        file<<le[j];
                        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;
// }
发表于 2010-2-1 13:36 | 显示全部楼层

回复 10楼 epistemer 的帖子

计算了好久,终于出结果了,贴出来大家看看对不对
数据: Lyapunov_Exponent.txt (24.12 KB, 下载次数: 18)

绘图:

绘图

绘图
发表于 2010-2-27 22:48 | 显示全部楼层
挺好的,算是基本算法了
发表于 2010-3-31 21:14 | 显示全部楼层

回复 12楼 xinwilliam 的帖子

怎么都是负的?
发表于 2010-4-1 08:49 | 显示全部楼层
图是怎么画的?matlab吗?
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-17 20:42 , Processed in 0.068939 second(s), 24 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表