wangjianbao11 发表于 2011-4-17 08:49

lorenz方程的求解问题

function ydot=hhfun(t,y,flag,h,a,b,r)
ydot=[-a*y(1)-a*y(2);
    -r*y(1)-y(1)-y(1)*y(3);
-b*y(3)+y(1)*y(2)];
%%program hh.m
h=0.0625;
a=10;b=10,r=40;
=ode23('hhfun',,,[],h,a,b,r);
figure
subplot(2,2,1);plot(t,y(:,1),t,y(:,2),t,y(:,3),'--')
subplot(2,2,2);plot(t,y(:,1),'b');Xlabel('t'),Ylabel('X');legend('x=5.001,y=z=5');
subplot(2,2,3);plot(t,y(:,2),'--');Xlabel('t'),Ylabel('Y');
subplot(2,2,4);plot(t,y(:,3),'--');Xlabel('t'),Ylabel('Z');

figure
subplot(2,2,1);plot(y(:,2),y(:,1));Xlabel('y与x');
subplot(2,2,2);plot(y(:,1),y(:,3));Xlabel('x与z');
subplot(2,2,3);plot(y(:,2),y(:,3));Xlabel('y与z');
subplot(2,2,4);plot(t,y(:,1),t,y(:,2),t,y(:,3),'--')
hold on
figure
axis ([-40 20 -40 80 -30 30])
plot3(y(:,1),y(:,2),y(:,3))

octopussheng 发表于 2011-4-17 18:59

这是？
问题还是分享？

meiyongyuandeze 发表于 2011-4-18 08:16

本帖最后由 meiyongyuandeze 于 2011-4-18 08:20 编辑

建议楼主把你的报错信息一给出，这样会省了很多时间。
刚才运行了下发现程序初始写的有问题，初值是列向量，应该是
=ode23('hhfun',,,[],h,a,b,r);

wangjianbao11 发表于 2011-4-26 15:14

我想问一下lorenz的分岔图，这是相图

meiyongyuandeze 发表于 2011-4-26 15:51

回复 4 # wangjianbao11 的帖子

Lorenz系统：
function dy = Lorenz(t,y)
% Lorenz系统
% 系统微分方程：
% dx/dt = -a(x-y)
% dy/dt = x(r-z)-y
% dz/dt = xy-bz
% a=y(4)
% r=y(5)
% b=y(6)
dy=zeros(6,1);
dy(1)=-y(4)*(y(1)-y(2));
dy(2)=y(1)*(y(5)-y(3))-y(2);
dy(3)=y(1)*y(2)-y(6)*y(3);
dy(4)=0;
dy(5)=0;
dy(6)=0;
function = getmax(y)
a=length(y);
j=1;
for i=(a-1)/2:a
b=(y(i,1)-y(i-2,1))/2;
c=(y(i,1)+y(i-2,1))/2-y(i-1,1);
if y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)&c==0
Xmax(j)=y(i-1,1);
j=j+1;
elseif y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)
Xmax(j)=y(i-1,1)-b^2/(4*c);
j=j+1;
end
end
function Lorenz_bifur_r_getmax
% 最大值法求解分岔图
clear all
t0=;%积分时间
%bifurcation
for r=linspace(1,500,1000); %r的变化精度
=ode45('Lorenz',t0,);
=getmax(y(:,1));
plot(r,Xmax,'b','markersize',1)
hold on
clear Xmax
end
给你发个论坛以前发过的一个用最大值方法计算分岔图的程序吧，你可以参照一下！

wangjianbao11 发表于 2011-4-26 16:19

谢谢您，师哥，运行结果和这个图形的分析你能给俺讲一下吗

meiyongyuandeze 发表于 2011-4-26 16:42

你先自己按照你的方程修改上面的程序，算出结果以后有了图才好分析！
上面程序的主要就是将不同参数时，稳态响应的最大值都找出来并画图，很好理解周期一就一个最大值，混沌有多个！

hsfy919 发表于 2011-4-27 13:09

回复 6 # wangjianbao11 的帖子

具体的非线性领域是相当复杂的，建议你先找一些书上的例子模仿一下

wangjianbao11 发表于 2011-4-27 13:43

unction dx=lorze1(t,X)
x=X(1);

y=X(2);

z=X(3);

a=10;

b=8/3;

global c;

dx=zeros(3,1);

dx(1)=-a*(x-y);

dx(2)=-x*z+c*x-y;

dx(3)=-b*z+x*y
clear;
global c;
range=;
k=0;
YY=[];
for c=range
    c
   y0=rand(1,3);
   k=k+1;
    tspan=;
   =ode45('lorze1',tspan,y0);
    count=find(t>100);
    Y=Y(count,:);
    % 画x的分岔图。
    j=1;
   n=length(Y(:,1));
    for i=2:n-1
      if Y(i-1,1)+eps<Y(i,1) && Y(i,1)>Y(i+1,1)+eps
% 简单的取出局部最大值。
            YY(k,j)=Y(i,1);
             j=j+1;
      end
   end
    if j>1
      plot(c,YY(k,),'k.','markersize',1);
   end
    hold on;
    index(k)=j-1;
end
xlabel('r');
ylabel('x max');

大家看看这程序有什么问题，运行出来结果我感觉不对

wangjianbao11 发表于 2011-4-27 13:48

C:\Users\fuyanzhong\Desktop\3.jpg

hsfy919 发表于 2011-4-27 14:37

回复 10 # wangjianbao11 的帖子

你上传的图片可能有问题，请重新编辑

wangjianbao11 发表于 2011-4-27 15:19

我也上传不上去，您能帮我做一个吗？谢谢了

wangjianbao11 发表于 2011-4-30 16:54



这是我用最大值法作出的图，其中a=10.b=8/3，c=, 各位大哥，给好好分析一下，看有什么错误
页: [1]
查看完整版本: lorenz方程的求解问题