MATLAB四阶龙格库塔方程求解方程组

Question

请问这个方程组如何用MATLAB的四阶龙格库塔方程求解？求一个代码，悬赏可加

Answer 1

（1）matlab基于四阶龙格库塔算法求解常微分方程组的函数ode45介绍
ode45是matlab中被广泛使用的基于四阶龙格库塔算法求解非线性常微分方程组的函数

同时，matlab中ode45的代码是开源的，可以直接在matlab主界面通过open ode45.m命令打开后拷贝：

（2）给出描述问题对应常微分方程组的函数文件
调用ode45时，需要先给出一个函数文件来描述常微分方程组，该函数的输入包括自变量t、当前各因变量的值、当前的参数值，输出是当前各因变量的导数值。相应的matlab代码如下（注意要新建一个与函数名一致的.m文件）：

function d_vars_div_dt=ode_fun(t,vars,b_x,g,b_y,k_z,Kzz)

%分配因变量的当前值
x=vars(1);y=vars(2);z=vars(3);v=vars(4);theta2=vars(5);
psi1=vars(6);phi1_prime=vars(7);phi_alpha_prime=vars(8);phi1=vars(9);phi_alpha=vars(10);
delta1=asin(sin(phi1)*cos(psi1)-sin(psi1)*cos(phi1)*cos(phi_alpha-theta2));
delta2=asin(sin(phi_alpha-theta2)*cos(psi1)/cos(delta1));

%计算当前个因变量的导数值（微商值）
dx_div_dt=v*cos(theta2)*cos(psi1);
dy_div_dt=v*sin(theta2)*cos(psi1);
dz_div_dt=v*sin(psi1);
dv_div_dt=-b_x*v^2-g*sin(theta2)*cos(psi1);
d_theta2_div_dt=(b_y*v^2*delta2-g*cos(theta2))/(v*cos(psi1));
d_psi1_div_dt=b_y*v*delta1+(g*sin(theta2)*sin(psi1))/v;
d_phi1_prime_div_dt=k_z*v^2*delta1-Kzz*v*phi1_prime;
d_phi_alpha_pirme_div_dt=k_z*v^2*delta2-Kzz*v*phi_alpha_prime;
d_phi1_div_dt=phi1_prime;
d_phi_alpha_div_dt=phi_alpha_prime;
d_vars_div_dt=real([dx_div_dt,dy_div_dt,dz_div_dt,dv_div_dt,d_theta2_div_dt,...
    d_psi1_div_dt,d_phi1_prime_div_dt,d_phi_alpha_pirme_div_dt,d_phi1_div_dt,d_phi_alpha_div_dt].');

end

（3）调用ode45求解问题中的常微分方程组

clear;clc;

%初始化参数和初值
b_x=1;g=2;b_y=3;k_z=4;Kzz=5;  %设置参数
t=(0:0.01:10).';  %设置时间采样点向量
init_x=rand();init_y=rand();init_z=rand();init_v=4;init_theta2=rand();
init_psi1=rand();init_phi1_prime=rand();init_phi_alpha_prime=rand();
init_phi1=rand();init_phi_alpha=rand();
init_vars=[init_x,init_y,init_z,init_v,init_theta2,...
    init_psi1,init_phi1_prime,init_phi_alpha_prime,init_phi1,init_phi_alpha];

%调用ode45函数求解非线性常微分方程组
[result_t,result_vars]=ode45(@(t,vars)ode_fun(t,vars,b_x,g,b_y,k_z,Kzz),t,init_vars);

%绘制求解结果
line_type_set={'r','g','b','k','m','r--','g--','b--','k--','m--'};
for i=1:size(result_vars,2)
    plot(result_t,result_vars(:,i),line_type_set{i},'linewidth',2);
    hold on;
end
axis tight;
xlabel('t');
ylabel('Amplitude');
h0=legend('x','y','z','v','\theta_2',...
    '\psi_1','\phi_1^''','\phi_{\alpha}^''','\phi_1','\phi_\alpha',...
    'location','bestoutside');
set(h0,'FontName','Times New Roman','FontSize',12,'FontWeight','bold','LineWidth',1);
h = gca;
set(h,'FontName','Times New Roman','FontSize',12,'FontWeight','bold','LineWidth',1);

（4）程序运行结果