Hello大家好!最近帮人家解决了一个问题,要求是这样的,要求用4阶龙格库塔法画分叉图。分叉混沌图在电力系统,动力系统的研究中比较广泛,常用来对系统的周期性进行一些研究。我在网上找了不少资料,没有实质性的进展,这方面的信息比较匮乏,于是我把我知道的解决方法做个总结,也算是学习记录,分享给大家。
其实四阶龙格库塔法只能得到一个迭代式,而分叉图,也是根据迭代式来画的。基本上有了迭代式就可以画出分叉混沌图了。 这里先给几个例子。
1、作分叉图的基本原理
1.1、第一个例子
%MATLAB代码
%要求:混沌与分叉 利用迭代格式x(k+1)=λsin(pi*x(k)),做出相应的Feigenbaum图。
clc,clear
y=@(k,x)k*sin(pi*x);
x0=0.3;
hold on
axis([1,3,-1,2]);
grid
for k=1:0.01:2;
for i=1:300
x0=y(k,x0);
if i>150
pause(0.01)
plot(k,x0,'.b')
% hold on;
end
end
end
grid
这个例子是用一个简单的迭代式来画分叉混沌图的,从这个例子里,不难学到如何如何画分叉图,其原理不难,两个循环来在2维平面作图,一个if筛选迭代的结果。
还是类似上面那个例子,这里给出不一样的代码实现方法。没看懂的可以再看一遍,里面主要的还是两个for循环,下面这个代码是用第三个for循环来实现对迭代结果的筛选的,这里不必纠结细节,功能实现用if或者for都可以。
clear;clf;
%feigenbaum曾对超越函数y=r*sin(pi*x)进行
%了分岔与混沌的研究,利用迭代
%格式x(k+1)=r*sin(pi*x)作出相应的feigenbaum图,并写出对应的matlab程序
hold on
axis([0,3,-3,3]);
grid
for a=0:0.005:3
x=[0.1];
for i=2:150
x(i)=a*sin(pi*x(i-1));
end
pause(0.1)
for i=101:150
plot(a,x(i),'k.');
end
end
1.2、第二个例子
这里再给另外一个方程的作图
%Henon mapping
%x(j+1)=1+yj-a*xj*xj
%y(j+1)=b*xj
clc
clear all
close all
a=0.4;b=0.3;
xjL=[];yjL=[];
for k=1:140
a=0.5+(k-1)*0.01;
xj=[];yj=[];
x=0.1;y=0.1;
for j=1:200
x=1+y-a*x*x;y=b*x;
xj=[xj;j x];
yj=[yj;j y];
if j>=120
xjL=[xjL;a b x];
yjL=[yjL;a b y];
end
end
end
plot(xjL(:,1),xjL(:,3),'k.')
xlabel('\itt','FontSize',20,'FontName','Times New Roman');
ylabel('\itx','FontSize',20,'FontName','Times New Roman');
grid on;
2、四阶龙格库塔法求微分方程
理论上,现在只要得到某个方程的迭代式似乎都可以能画出它的分叉图,而4阶龙格库塔法只需要套一下公式就好了:
%%目标高阶常微分线性方程:mx''+c*x'+k*x+u*x*x+j*x^3 = F*cos(w*t)
for c=0:0.003:1;
fun=@(c,x)[x(2);(-c.*x(2)-k.*x(1)-u.*x(1).^2-j.*x(1).^3+F.*cos(w.*t))./m];
for n = 1:200
k1 = fun(c,x(n,:)');
k2 = fun(c+h/2,x(n,:)'+h/2*k1);
k3 = fun(c+h/2,x(n,:)'+h/2*k2);
k4 = fun(c+h,x(n,:)'+h*k3);
x(n+1,:) = x(n,:)+h/6*(k1+2*k2+2*k3+k4)';
end
end
这里只是我的方程解法,网上这方面的例子有很多。我这里再给一个以前上课时老师给的例子吧:
function R=rk4(f,a,b,ya,N)
%y'=f(x,y)
%a,b左右端点。
%N为迭代步数。
%h为步长。
%ya为初值。
h=(b-a)/N;
T=zeros(1,N+1);
Y=zeros(1,N+1);
T=a:h:b;
Y(1)=ya;
for j=1:N
k1=h*feval(f,T(j),Y(j));
k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);
k3=h*feval(f,T(j)+h/2,Y(j)+k2/2);
k4=h*feval(f,T(j)+h,Y(j)+k3);
Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;
end
R=[T' Y']
T
Y
另外啰嗦一句,MATLAB常用的解微分方程有ode,欧拉方法,龙格库塔,Dsolve。
3、最后我的代码与结果
MATLAB代码
clc;
clear all;
close all;
%% mx''+c*x'+k*x+u*x*x+j*x^3 = F*cos(w*t),
m=1;
k=90;%80~110之间出现分叉,越大图越偏左移
u=-15;%正负数绝对值数值趋势一样,正数范围内,数字越小图越左越大,负数时,数越小图越小密集图越偏右
j=60;%正数数值越大,图越大,-12图消失,70四分叉消失,40到70之间
F=115;
w=pi;
%%%%%%%
t0=2;
%tf=300;
h=0.1;
t=t0;
%%%%%%%%%%
%c=0:0.01:1;
%%%%%%%
x(1,1)=t;%t初值对数值很敏感,直接关系到分叉的出现
x(1,2)=6;%与分支长短有关,对数据敏感,敏感程度略低于t0
%%mx''+c*x'+k*x+u*x*x+j*x^3 = F*cos(w*t),0<=t<=1
%for t=0.05:0.01:1%?á×?±ê
hold on
axis([0,1,-1.1,2]);
grid
for c=0:0.003:1;%1分4在0.001时边缘可见
fun=@(c,x)[x(2);(-c.*x(2)-k.*x(1)-u.*x(1).^2-j.*x(1).^3+F.*cos(w.*t))./m];
for n = 1:200 %%用4阶龙格库塔法求解微分方程得到迭代式,有了迭代式就可以作分叉混沌图了
k1 = fun(c,x(n,:)');
k2 = fun(c+h/2,x(n,:)'+h/2*k1); % ????????±?×????????ò???à??
k3 = fun(c+h/2,x(n,:)'+h/2*k2);
k4 = fun(c+h,x(n,:)'+h*k3);
x(n+1,:) = x(n,:)+h/6*(k1+2*k2+2*k3+k4)';%???????×Runge-Kutta????++
% plot(c(n),k1(1,1),'-k.')
if n>100 %%单一个横坐标时对纵坐标的值进行筛选,也就是迭代次数大于100的才可以画在图上,迭代次数越多,精度越高,点也会越少
pause(0.00001) %这一步是为了看到作图过程,类似短暂的暂停画布时间
% cc=1-c;
plot(c,x(n+1),'.b') %%所谓的分叉混沌图,本质上是两个循环,外循环是横坐标,内循环是纵坐标的值,这里的内循环的值通过4阶龙格库塔法迭代得到
%plot(c,k1(1,1),'.k')
end
% hold on
end
end