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:300x0=y(k,x0);if i>150pause(0.01)plot(k,x0,'.b')% hold on;endend 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:3x=[0.1]; for i=2:150x(i)=a*sin(pi*x(i-1)); endpause(0.1)for i=101:150plot(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:140a=0.5+(k-1)*0.01;xj=[];yj=[];x=0.1;y=0.1;for j=1:200x=1+y-a*x*x;y=b*x;xj=[xj;j x]; yj=[yj;j y];if j>=120xjL=[xjL;a b x]; yjL=[yjL;a b y]; endend 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 onend end

總結

以上是生活随笔為你收集整理的利用四阶龙格库塔法（Runge-Kutta methods）求解常微分方程并用其迭代式用MATLAB绘制分叉混沌图的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。