¸÷λÂÛ̳µÄ´óÉñ£¬Ð¡µÜ±¾¿Æ±ÏÒµÂÛÎÄ£¬ÀÏʦ¸øµÄÄ£ÐÍÊÇÕâÑùÒ»¸öº¬¼ä϶µÄ¶þ×ÔÓɶÈÅöײϵͳ£¬Ñо¿ËüµÄ·Ö²æÓë»ìã磺×óÓÒ¸÷Á½¸öСÇò£¬¸÷×ÔÁ¬½Ó×ÅÒ»¸ö·ÇÏßÐԵĵ¯»ÉºÍ×èÄáÆ÷£¬¼ÓÉϼ¤ÀøÁ¦
¾ßÌåµÄÎÒÉÏ´«ÁËͼƬ£¬»¹Óи½¼þ¡£
ÎÊÌâÃèÊöÍê±Ï£¬ÎÒÏÖÔÚÓöµ½µÄÎÊÌâÊÇ£ºÒªÔÚmatlabÀÂÇÕâÁ½¸öСÇòµÄÅöײ£¬ÀÏʦ˵Ҫ°Ñʱ¼ä²½ÉèÖÃΪ1e-6pi£¬¼ÆËãʱ¼äΪ500pi£¬²¢ÇÒÈ¥µôÇ°Ãæ200piµÄʱ¼ä¡£µ«ÊÇ£¡Õâ¸ö¸ù±¾Ëã²»ÁËÍÛ£¬ÎÒÕâ¸öµçÄÔËãÕâô¶àµÄ»°µÃËãµ½ºïÄêÂíÔÂÈ¥ÁË£¬È»ºóÎÒ×ԾͰÑʱ¼äËõ¶Ìµ½50pi£¬Ê±¼ä²½¸ÄΪ1e-3pi£¬µ«ÊǼÆËãÆðÀ´Ò²·Ñʱ¼ä£¬ÎÒÏÖÔÚ¾ÍÏëÎÊÎʸ÷λ£º
1.Õâ¸öÎÊÌâ±ØÐëÒªÕâô³¤µÄʱ¼äÂð£¿ÓÐʲô°ì·¨¿ÉÒÔËõ¼õһϣ¿
2.ÔõÑùÅжÏÕâÁ½¸öСÇòÅöײÉÏÁË£¿ ÎÒдµÄ´úÂëÊÇ£ºif abs(delta-abs(x(1,n+1)-x(2,n+1)))<=10^(-2) deltaÊǼä϶µÄ¾àÀë¡£x£¨1£¬n+1£©ÊǵÚn+1²½Ê±µÚÒ»¸öÇòµÄλÖã¬x£¨2£¬n+1£©ÊǵÚn+1²½Ê±µÚ¶þ¸öÇòµÄλÖᣲ»ÖªµÀÕâÑùд¶Ô²»¶Ô¡£
×òÍíËãµÄ½á¹ûÎÒ½ØÁËͼ£¬´«¸½¼þÁË¡£ Ï£Íû´ó¸çÃÇÄܰïСµÜ¿´¿´¡£
¶ÔÁË£¬»¹ÓÐÎÒ×Ô¼º±àµÄmatlab´úÂ룺
µÚÒ»¸ödyfunº¯Êý£¬ÓÃÓÚ´æ·Å¿ØÖÆ·½³Ì
function f=dyfun(x,y)
f(1)=y(3);
f(2)=y(4);
f(3)=-0.01*y(3)*0.1-y(1)*1^2-0.01*y(1)^3+10*0.01*cos(x*2);
f(4)=-0.01*y(4)*0.15-y(2)*1.5^2-0.012*y(2)^3+3*0.01*cos(x*2);
f=f( : );
end
--------------------------------------------------------------------------------------------------
µÚ¶þ¸öº¯Êý£¬×Ô¼º±àдµÄ4½×R-K·½·¨£¬ÒòΪode45ÀïÃæÎÒ²»ÖªµÀÔõôÌí¼ÓÅжÏÌõ¼þ£¬ÅжÏÁ½¸öСÇòÊÇ·ñײÉÏÁË¡£
function[t,x]=narunge4s(dyfun,tspan,x0,h)
delta=input('delta=');
t=tspan(1):h:tspan(2);
x=zeros(length(x0),length(t));
x(:,1)=x0(: );
for n=1: (length(t)-1)
k1=feval(dyfun,t(n),x(:,n));
k2=feval(dyfun,t(n)+h/2,x(:,n)+h/2*k1);
k3=feval(dyfun,t(n)+h/2,x(:,n)+h/2*k2);
k4=feval(dyfun,t(n+1),x(:,n)+h*k3);
x(:,n+1)=x(:,n)+(h/6)*(k1+2*k2+2*k3+k4);
%ÏÂÃæÕâÈýÐÐÊÇÓÉÅöײ¹ý³ÌÖж¯Á¿ÊغãºÍ»Ö¸´ÏµÊýµÃµ½µÄ
% lanbda=n1/n2 ´Ë´¦È¡×ö1
% q+p=x(3,n+1)+x(4,n+1);
% q-p=-0.8*(x(3,n+1)-x(4,n+1));
if abs(delta-abs(x(1,n+1)-x(2,n+1)))<=10^(-2)
%Õâ¾ÍÊÇÅжÏÊÇ·ñÅöײµÄÓï¾ä£¬²»ÖªµÀ¶Ô²»¶Ô
q=(x(3,n+1)+x(4,n+1)-0.8*(x(3,n+1)-x(4,n+1)))/2;
p=(x(3,n+1)+x(4,n+1)+0.8*(x(3,n+1)+x(4,n+1)))/2;
x(3,n+1)=q;
x(4,n+1)=p; Õâ¾ÍÊÇÅжÏÊÇ·ñÅöײµÄÓï¾ä£¬²»ÖªµÀ¶Ô²»¶Ô
% stop=input('stop=');
end
n
% ÕâÊÇΪÁËÒ»±ß¼ÆËãÒ»±ß¿´Í¼£¬¸ô100¸öµã»Ò»´Î
if mod(n,100)==0
figure(1);
plot(t(1:n),x(1,1:n),'r',t(1:n),x(2,1:n),'r-.');
legend('x(1)','x(2)');
hold on
xlabel('x(1)');
ylabel('x(2)');
figure(2);
plot(x(1,1:n),x(3,1:n),'r');
xlabel('x(1)');
ylabel('dx(1)/dt');
title('The 1st')
figure(3);
plot(x(2,1:n),x(4,1:n),'r-')
xlabel('x(2)');
ylabel('dx(2)/dt');
title('The 2nd')
hold on
end
end
-------------------------------------------------------------------------------------------------
Ö÷³ÌÐò£¬µ÷ÓÃ֮ǰµÄ¼ÆËã¡£
clc
clear all
[t,x]=marunge4s(@dyfun,[0 50*pi],[0 0.5 0 0],10^-4*pi);
![Çó°ïÖú£¡£¡¹ØÓÚÀûÓÃmatlabÄ£Äâ¶þ×ÔÓɶȺ¬¼ä϶ÅöײϵͳµÄÎÊÌ⣬±¾¿Æ±ÏÒµÂÛÎÄ£¡£¡£¡£¡£¡]()
ÂÛÎÄÄ£ÐÍ.PNG
![Çó°ïÖú£¡£¡¹ØÓÚÀûÓÃmatlabÄ£Äâ¶þ×ÔÓɶȺ¬¼ä϶ÅöײϵͳµÄÎÊÌ⣬±¾¿Æ±ÏÒµÂÛÎÄ£¡£¡£¡£¡£¡-1]()
Á½¸öСÇòµÄÎ»ÒÆ-ʱ¼äͼ.PNG
![Çó°ïÖú£¡£¡¹ØÓÚÀûÓÃmatlabÄ£Äâ¶þ×ÔÓɶȺ¬¼ä϶ÅöײϵͳµÄÎÊÌ⣬±¾¿Æ±ÏÒµÂÛÎÄ£¡£¡£¡£¡£¡-2]()
µÚ¶þ¸öСÇòµÄÎ»ÒÆ-ËÙ¶ÈÏàͼ.PNG |