请问用四阶龙格库塔法解二阶微分方程的思想是什么?
请问用四阶龙格库塔法解二阶微分方程的思想是什么?
最近我遇到了一个难题!就是求一个二阶微分方程,形式如:y''=f(y),初始条件是y(0)=0,y'(0)=0,y是t的函数.要求用四阶龙格库塔法求解!请问思路是什么?由于f(y)函数复杂,所以解析解是求不出来的,必须用四阶龙格库塔法解之,如果您有实例,希望能够和我分享一下!
令y1=y';y2=y''=(y1)';得到了二阶微分方程的一阶形式,然后可以按照普通的一节微分方程组的四阶龙格库塔法求解了;这样说的有点抽象,你可以把方程复制出来 我有空帮你编写一下好的 非常感谢呢
292200y''=k*(90-y)-8390.7*sin(y);y(0)=0,y'(0)=0;
其中y是时间t的函数,k是自己可以调的一个系数;步长取0.001吧,最少计算一千次,因为这是折叠机翼的转动微分方程,要求在1秒内展开,所以最少计算一千次。解这个微分方程的目的是选取k,从而实现1秒内展开,即y(t=1)>=90度
非常感谢哈
默认y的单位是弧度
k=1000;
t=0:0.001:1;
Y=[];
err=1
K=[];
Ymax=[];
xishu=1.01;
while err
X=[0 0];
k=xishu*k;
K=[K;k];
Y=[];
for i=1:1001
Y_1 = Runge_Kutta41(t(i),X,@folded_wing,0.001,k);
Y=[Y;t(i),Y_1];
X=Y_1;
end
ymax=max(Y(:,2));
Ymax=[Ymax;ymax];
if ymax>=pi/2&&ymax<=92/rad2deg(1)
break
elseif ymax>92/rad2deg(1)
xishu=0.98;
else
xishu=1.01;
end
end
plot(Y(:,1),Y(:,2),'-b','linewidth',2);
hold on
grid on
xlabel('t-time','fontsize',14)
ylabel('y(rad)','fontsize',14)
plot([0,1],[pi/2 pi/2],'-g','linewidth',2);
tstring=['292200y"=',num2str(k),'*(90-y)-8390.7*sin(y);'];
title(tstring,'fontsize',14);
legend('y"','y=pi/2')
function Y = Runge_Kutta41(t,X,f,h,k)
K1=f(t,X,k);
K2=f(t+h/2,X+h/2*K1,k);
K3=f(t+h/2,X+h/2*K2,k);
K4=f(t+h,X+h*K3,k);
Y=X+h/6*(K1+2*K2+2*K3+K4);
function dy=folded_wing(t,y,k)
% y1=y;y2=y'
dy=zeros(1,2);
dy(1)=y(2);
dy(2)=1/292200*(k*(90-y(1))-8390.7*sin(y(1)));