这种图能说明方程是刚性的吗？

zhailiangjun

这是我最近计算所得到的一些图像，是一个4自由度微分方程组的4个变量随时间变化的曲线图，很有意思。请大家看一下吧，这些能说明这个方程组是刚性的吗？

无水1324

确实 很奇怪的 ，你 看一下 方程 中的 x（3），x（4）是否写错

octopussheng

这个曲线确实很有意思，关注一下！

咕噜噜

有可能是你程序有问题，以前我得程序写错的时候就出现过此类情况:loveliness:

zhailiangjun

嗯，你说的错误是指那些方面的错误呐，我把我的程序贴上来，请大家看一下啊，:loveliness:
谢谢了。
clear all;
clc;
syms x4;
A1=-18.2219; A12=-2.85; N=42; beta=2.354; phi0=1.828; lamta=1.0571; u1=1.0097;
u3=1/15.9994; r1e=0.9558;     %一些参数

x1= (1.0-0.9)*rand(1,1)+0.9;
x2=(400-200)*rand(1,1)+200;
x3=(1.0-0.9)*rand(1,1)+0.9;

h=vpa((A1+A12)*N^2*(2-exp(-beta*(x1-r1e)))*exp(-beta*(x1-r1e))+(A1+A12)*N^2*(2-exp(-beta*(x3-r1e)))*exp(-beta*(x3-r1e))+2*A12*N*N*((2-exp(-beta*(x1-r1e)))*exp(-beta*(x1-r1e))*(2-exp(-beta*(x3-r1e)))*exp(-beta*(x3-r1e)))^(1/2)+1/4*lamta*N*N*(2*exp(-beta*(x1-r1e))+2*exp(-beta*(x3-r1e))-2*exp(-beta*(x1-r1e)-beta*(x3-r1e))-2*((2-exp(-beta*(x1-r1e)))*exp(-beta*(x1-r1e))*(2-exp(-beta*(x3-r1e)))*exp(-beta*(x3-r1e)))^(1/2))+1/2*((u1+u3)*x2^2+(u1+u3)*x4^2+2*u3*(cos(phi0)*x2*x4))-10287.2);
%这是一个代数哈密顿方程。
s=solve(h,'x4');
d=double(s)
e=d(1,1)             %随即选取初始点
%==============================================================
x11=x1;
x22=x2;
x33=x3;
x44=e;
[t,x]=ode113(@gudingfun,[0 8],[x11 x22 x33 x44]);     %计算微分方程
%================================================================
figure(1)
plot(t(:),x(:,1))
figure(2)
plot(t(:),x(:,2))
figure(3)
plot(t(:),x(:,3))
figure(4)
plot(t(:),x(:,4))

[ 本帖最后由 zhailiangjun 于 2008-4-11 10:40 编辑 ]

zhailiangjun

呵呵，谢谢你的关注啊，我把程序贴出来了，有空帮忙看一下啊。谢谢了。:lol

zhailiangjun

嗯，我仔细的检查一下，谢谢你的意见。呵呵。我把我的程序贴了上来，有空帮忙看一下哈，谢谢你了。无水大哥。

zhailiangjun

无水大哥，我现在要在进行贴图，该怎么做啊？:@L

无水1324

“发表新回复”在帖子的又下角，进去之后可以加附加的回复了

咕噜噜

子程序那？

zhailiangjun

谢谢你的回复，现在把子程序贴在这。有空帮帮忙啊，谢谢。
function dx=gudingfun(t,x);
u1=1/1.0081;
u3=1/16.0014;
A1=-18.2219;
A12=-2.85;
N=42;
beta=2.354;
phi0=1.828;
lamta=1.0571;
r1e=0.9558;
dx=zeros(4,1);

dx(1)=(u1+u3)*x(2)+u3*cos(phi0)*x(4);

dx(2)=-(A1+A12)*N^2*beta*exp(-beta*(x(1)-r1e))^2+(A1+A12)*N^2*(2-exp(-beta*(x(1)-r1e)))*beta*exp(-beta*(x(1)-r1e))-A12*N^2/((2-exp(-beta*(x(1)-r1e)))*exp(-beta*(x(1)-r1e))*(2-exp(-beta*(x(3)-r1e)))*exp(-beta*(x(3)-r1e)))^(1/2)*(beta*exp(-beta*(x(1)-r1e))^2*(2-exp(-beta*(x(3)-r1e)))*exp(-beta*(x(3)-r1e))-(2-exp(-beta*(x(1)-r1e)))*beta*exp(-beta*(x(1)-r1e))*(2-exp(-beta*(x(3)-r1e)))*exp(-beta*(x(3)-r1e)))-1/4*lamta*N^2*(-2*beta*exp(-beta*(x(1)-r1e))+2*beta*exp(-beta*(x(1)-r1e)-beta*(x(3)-r1e))-1/((2-exp(-beta*(x(1)-r1e)))*exp(-beta*(x(1)-r1e))*(2-exp(-beta*(x(3)-r1e)))*exp(-beta*(x(3)-r1e)))^(1/2)*(beta*exp(-beta*(x(1)-r1e))^2*(2-exp(-beta*(x(3)-r1e)))*exp(-beta*(x(3)-r1e))-(2-exp(-beta*(x(1)-r1e)))*beta*exp(-beta*(x(1)-r1e))*(2-exp(-beta*(x(3)-r1e)))*exp(-beta*(x(3)-r1e))));

dx(3)=(u1+u3)*x(4)+u3*cos(phi0)*x(2);

dx(4)=-(A1+A12)*N^2*beta*exp(-2*beta*(x(3)-r1e))+(A1+A12)*N^2*(2-exp(-beta*(x(3)-r1e)))*beta*exp(-beta*(x(3)-r1e))-A12*N^2/((2-exp(-beta*(x(1)-r1e)))*(2-exp(-beta*(x(3)-r1e)))*exp(-beta*(x(1)-r1e)-beta*(x(3)-r1e)))^(1/2)*((2-exp(-beta*(x(1)-r1e)))*beta*exp(-beta*(x(1)-r1e)-2*beta*(x(3)-r1e))+(-2+exp(-beta*(x(3)-r1e)))*beta*(2-exp(-beta*(x(1)-r1e)))*exp(-beta*(x(1)-r1e)-beta*(x(3)-r1e)))-1/4*lamta*N^2*(-2*beta*exp(-beta*(x(3)-r1e))+2*beta*exp(-beta*(x(1)-r1e)-beta*(x(3)-r1e))-1/((2-exp(-beta*(x(1)-r1e)))*(2-exp(-beta*(x(3)-r1e)))*exp(-beta*(x(1)-r1e)-beta*(x(3)-r1e)))^(1/2)*((2-exp(-beta*(x(1)-r1e)))*beta*exp(-beta*(x(1)-r1e)-2*beta*(x(3)-r1e))+(-2+exp(-beta*(x(3)-r1e)))*beta*(2-exp(-beta*(x(1)-r1e)))*exp(-beta*(x(1)-r1e)-beta*(x(3)-r1e))));

zhailiangjun

还是这个问题，前两天忽然又得到了这么一组图，请大家帮忙分析一下原因哈，
这个时候4个数据均出现了震荡解，但是只是在那么一组初始点数据里面比较合适。不知道为什么？:@L :@Q

zhailiangjun

这个就是我在计算的微分方程组了，呵呵。:@P 贴出来一个简单的形式请大家帮忙看一下这里面有什么问题吧。:lol

[ 本帖最后由 无水1324 于 2008-4-14 18:40 编辑 ]

无水1324

没有看出什么问题，但是你上面那个图应该是混沌的了

octopussheng

这个是时间历程曲线？？