lorenz方程的求解问题

wangjianbao11

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;
[t,y]=ode23('hhfun',[0:0.1:20],[5 5 5],[],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

这是？
问题还是分享？

meiyongyuandeze

建议楼主把你的报错信息一给出，这样会省了很多时间。
刚才运行了下发现程序初始写的有问题，初值是列向量，应该是

1. [t,y]=ode23('hhfun',[0:0.1:20],[5;5;5],[],h,a,b,r);

复制代码

wangjianbao11

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

meiyongyuandeze

回复 4 # wangjianbao11 的帖子

1. Lorenz系统：
   
2. function dy = Lorenz(t,y)
   
3. % Lorenz系统
   
4. % 系统微分方程：
   
5. % dx/dt = -a(x-y)
   
6. % dy/dt = x(r-z)-y
   
7. % dz/dt = xy-bz
   
8. % a=y(4)
   
9. % r=y(5)
   
10. % b=y(6)
    
11. dy=zeros(6,1);
    
12. dy(1)=-y(4)*(y(1)-y(2));
    
13. dy(2)=y(1)*(y(5)-y(3))-y(2);
    
14. dy(3)=y(1)*y(2)-y(6)*y(3);
    
15. dy(4)=0;
    
16. dy(5)=0;
    
17. dy(6)=0;
    
18. function [Xmax] = getmax(y)
    
19. a=length(y);
    
20. j=1;
    
21. for i=(a-1)/2:a
    
22. b=(y(i,1)-y(i-2,1))/2;
    
23. c=(y(i,1)+y(i-2,1))/2-y(i-1,1);
    
24. if y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)&c==0
    
25. Xmax(j)=y(i-1,1);
    
26. j=j+1;
    
27. elseif y(i-2,1)<=y(i-1,1)&y(i-1,1)>=y(i,1)
    
28. Xmax(j)=y(i-1,1)-b^2/(4*c);
    
29. j=j+1;
    
30. end
    
31. end
    
32. function Lorenz_bifur_r_getmax
    
33. % 最大值法求解分岔图
    
34. clear all
    
35. t0=[0 100];%积分时间
    
36. %bifurcation
    
37. for r=linspace(1,500,1000); %r的变化精度
    
38. [t,y]=ode45('Lorenz',t0,[1;1;1;16;r;4]);
    
39. [Xmax]=getmax(y(:,1));
    
40. plot(r,Xmax,'b','markersize',1)
    
41. hold on
    
42. clear Xmax
    
43. end
    

复制代码

给你发个论坛以前发过的一个用最大值方法计算分岔图的程序吧，你可以参照一下！

wangjianbao11

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

meiyongyuandeze

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

hsfy919

回复 6 # wangjianbao11 的帖子

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

wangjianbao11

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=[0:0.1:16];
k=0;
YY=[];
for c=range
    c
     y0=rand(1,3);
     k=k+1;
    tspan=[0:0.1:200];
     [t,Y]=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,[1:j-1]),'k.','markersize',1);
     end
    hold on;
    index(k)=j-1;
end
xlabel('r');
ylabel('x max');

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

wangjianbao11

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

hsfy919

回复 10 # wangjianbao11 的帖子

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

wangjianbao11

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

wangjianbao11

这是我用最大值法作出的图，其中a=10.b=8/3，c=[0,30], 各位大哥，给好好分析一下，看有什么错误