这种图能说明方程是刚性的吗？
这是我最近计算所得到的一些图像，是一个4自由度微分方程组的4个变量随时间变化的曲线图，很有意思。请大家看一下吧，这些能说明这个方程组是刚性的吗？回复 楼主 的帖子
确实 很奇怪的 ，你 看一下 方程 中的 x（3），x（4）是否写错 这个曲线确实很有意思，关注一下！ 有可能是你程序有问题，以前我得程序写错的时候就出现过此类情况:loveliness:回复 4楼 的帖子
嗯，你说的错误是指那些方面的错误呐，我把我的程序贴上来，请大家看一下啊，: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;
=ode113(@gudingfun,,); %计算微分方程
%================================================================
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 编辑 ]
回复 3楼 的帖子
呵呵，谢谢你的关注啊，我把程序贴出来了，有空帮忙看一下啊。谢谢了。:lol回复 2楼 的帖子
嗯，我仔细的检查一下，谢谢你的意见。呵呵。我把我的程序贴了上来，有空帮忙看一下哈，谢谢你了。无水大哥。回复 2楼 的帖子
无水大哥，我现在要在进行贴图，该怎么做啊？:@L回复 8楼 的帖子
“发表新回复”在帖子的又下角，进去之后可以加附加的回复了 子程序那？回复 10楼 的帖子
谢谢你的回复，现在把子程序贴在这。有空帮帮忙啊，谢谢。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))));
新的结果图
还是这个问题，前两天忽然又得到了这么一组图，请大家帮忙分析一下原因哈，这个时候4个数据均出现了震荡解，但是只是在那么一组初始点数据里面比较合适。不知道为什么？:@L :@Q
微分方程组
这个就是我在计算的微分方程组了，呵呵。:@P 贴出来一个简单的形式请大家帮忙看一下这里面有什么问题吧。:lol
[ 本帖最后由 无水1324 于 2008-4-14 18:40 编辑 ]