声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1706|回复: 1

[计算数学] 求大神把程序耦合一下,或者看看问题出在哪

[复制链接]
发表于 2013-9-12 11:33 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
下面是计算分数阶吸引子的程序,分数阶的系统已经计算出来了,按照整数阶的规律画李指数应该很简单,可惜不会。顺便附带整数阶李指数的程序,看有木有人能给耦合一下,
    function [T, Y]=FOChen(parameters, orders, TSim, Y0)

% time step:
h=0.005;
% number of calculated mesh points:

n=round(TSim/h);

%orders of derivatives, respectively:

q1=orders(1); q2=orders(2); q3=orders(3);

% constants of Chen's system:
r=parameters(1);
% binomial coefficients calculation:
cp1=1; cp2=1; cp3=1;
for j=1:n
     c1(j)=(1-(1+q1)/j)*cp1;
     c2(j)=(1-(1+q2)/j)*cp2;
     c3(j)=(1-(1+q3)/j)*cp3;
     cp1=c1(j); cp2=c2(j); cp3=c3(j);
end
% initial conditions setting:
  x(1)=Y0(1); y(1)=Y0(2); z(1)=Y0(3);
% calculation of phase portraits /numerical solution/:
for i=2:n
     x(i)=(10*(y(i-1)-x(i-1)))*h^q1 - memo(x, c1, i);
     y(i)=(-x(i)*z(i-1)+(24-r*4)*x(i)+r*y(i-1))*h^q2 - memo(y, c2, i);
     z(i)=(x(i)*y(i)-8/3*z(i-1))*h^q3 - memo(z, c3, i);
end
for j=1:n
     Y(j,1)=x(j);
     Y(j,2)=y(j);
     Y(j,3)=z(j);
end
T=h:h:TSim;


[t,y]=FOchen([5], [0.95 0.95 0.95], 100, [-8.3 -10  12]);
figure
plot3(y(:,1), y(:,3), y(:,2), 'k');
xlabel('z(t)'); ylabel('x(t)'); zlabel('y(t)'); grid;
title('chenxiyinzi x,y,z')



Z1=[];Z2=[];Z3=[];
global a;
global b;
global r;
a=10;b=8/3;
for r=linspace(0,5,40);
    y=[1;1;1;1;0;0;0;1;0;0;0;1];
    lp=0;
    for k=1:100
      [T,Y]=FOchen([5], [0.95 0.95 0.95], 10, [-8.3 -10  12]);
        y = Y(size(Y,1),:);
        y0 = [y(4) y(7) y(10);
            y(5) y(8) y(11);
            y(6) y(9) y(12)];
        y0=GS(y0);
        mod(1)=norm(y0(:,1));
        mod(2)=norm(y0(:,2));
        mod(3)=norm(y0(:,3));
        lp = lp+log(abs(mod));
        y0(:,1)=y0(:,1)/mod(1);
        y0(:,2)=y0(:,2)/mod(2);
        y0(:,3)=y0(:,3)/mod(3);
        y(4:12) = y0';
    end
    lp=lp/100;
    Z1=[Z1 lp(1)];
    Z2=[Z2 lp(2)];
    Z3=[Z3 lp(3)];
end
r=linspace(0,5,40);
plot(r,Z1,'-',r,Z2,'-',r,Z3,'-');
title('Lyapunov exponents of Lorenz')
xlabel('parameter r'),ylabel('lyapunov exponents')
grid on

正交化;
function A = GS(V)  
v1 = V(:,1);
v2 = V(:,2);
v3 = V(:,3);
a1 = zeros(3,1);
a2 = zeros(3,1);
a3 = zeros(3,1);
a1 = v1;
a2 = v2-((a1'*v2)/(a1'*a1))*a1;
a3 = v3-((a1'*v3)/(a1'*a1))*a1-((a2'*v3)/(a2'*a2))*a2;
A = [a1,a2,a3];

回复
分享到:

使用道具 举报

 楼主| 发表于 2013-9-20 14:34 | 显示全部楼层
     大神们,想求个做分数阶的李指数,开价300人人民币,,大家讨论下吧,
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-10 17:47 , Processed in 0.063689 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表