|
|
这段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;
// } |
|
发表于 2006-6-17 18:13
提升卡
变色卡
显身卡
发表于 2006-7-13 22:04
