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

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

（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）程序运行结果

可以编写四个子函数如下：

function output = f1(x,u1,u2,w1,w2)
output = u2;

function output = f2(x,u1,u2,w1,w2)
output = cos(x)-u2+w2;

function output = f3(x,u1,u2,w1,w2)
output = w2;

function output = f4(x,u1,u2,w1,w2)
output = -sin(x)-2*cos(x)-w2+u2;

四阶龙格库塔函数程序

function [u1,u2,w1,w2] = RK4_2variable(u1,u2,w1,w2,h,a,b)

x = a:h:b;

for i = 1:length(x)-1

k11 = f1(x(i) , u1(i) , u2(i) , w1(i) , w2(i));
k21 = f2(x(i) , u1(i) , u2(i) , w1(i) , w2(i));
L11 = f3(x(i) , u1(i) , u2(i) , w1(i) , w2(i));
L21 = f4(x(i) , u1(i) , u2(i) , w1(i) , w2(i));

k12 = f1(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2);
k22 = f2(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2);
L12 = f3(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2);
L22 = f4(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2);

k13 = f1(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2);
k23 = f2(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2);
L13 = f3(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2);
L23 = f4(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2);

k14 = f1(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);
k24 = f2(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);
L14 = f3(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);
L24 = f4(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);

u1(i+1) = u1(i) + h/6 * (k11 + 2*k12 + 2*k13 + k14);
u2(i+1) = u2(i) + h/6 * (k21 + 2*k22 + 2*k23 + k24);
w1(i+1) = w1(i) + h/6 * (L11 + 2*L12 + 2*L13 + L14);
w2(i+1) = w2(i) + h/6 * (L21 + 2*L22 + 2*L23 + L24);

end
end

计算主程序：

% main func
clear;clc;
u1(1) = 0;
u2(1) = 1;
w1(1) = 1;
w2(1) = 0;
h=0.01;
a = 0;b=20;
[u1,u2,w1,w2] = RK4_2variable(u1,u2,w1,w2,h,a,b);
figure 
plot(a:h:b,u1,'r-');
hold on
plot(a:h:b,sin(a:h:b),'b-.');
xlabel('time');
ylabel('x');
legend('计算值','精确值');

（望采纳）
1.四阶龙格库塔函数程序，变量命名可以对照公式看一下

function [xout, yout, zout, vout, theta2out, bfai1out] = Runge_Kutta4(g,bx,by,x0,y0,z0,v0,theta20,bfai10,sfai1,faiapha,h,ta,tb,n)
% ta和tb是时间范围的左值和右值，n是时间的微分间隔 ，其余的就是参数和初值了
xout = [x0];
yout = [y0];
zout = [z0];
vout = [v0];
theta2out = [theta20];
bfai1out = [bfai10];

xn = x0;
yn = y0;
zn = z0;
vn = v0;
theta2n = theta20;
bfai1n = bfai10;


for tn = ta:n:tb
zigama1 = asin(sin(sfai1)*cos(bfai1n)-sin(bfai1n)*cos(sfai1)*cos(faiapha-theta2n));
zigama2 = asin(sin(faiapha-theta2n)*cos(bfai1n)/cos(zigama1));
k1_x = h*vn*cos(theta2n)*cos(bfai1n);
k1_y = h*vn*sin(theta2n)*cos(bfai1n);
k1_z = h*vn*sin(bfai1n);
k1_v = h*(-bx)*vn^2-g*sin(theta2n)*cos(bfai1n);
k1_theta2 = h*((by*vn^2*zigama2-g*cos(theta2n))/(vn*cos(bfai1n)));
k1_bfai1 = h*(by*vn*zigama1 + (g*sin(theta2n)*sin(bfai1n)/vn));

k2_x = h*(vn+k1_v/2)*cos(theta2n+k1_theta2/2)*cos(bfai1n+k1_bfai1/2);
k2_y = h*(vn+k1_v/2)*sin(theta2n+k1_theta2/2)*cos(bfai1n+k1_bfai1/2);
k2_z = h*(vn+k1_v/2)*sin(bfai1n+k1_bfai1/2);
k2_v = h*(-bx)*(vn+k1_v/2)^2-g*sin(theta2n+k1_theta2/2)*cos(bfai1n+k1_bfai1/2);
k2_theta2 = h*((by*(vn+k1_v/2)^2*zigama2-g*cos(theta2n+k1_theta2/2))/((vn+k1_v/2)*cos(bfai1n+k1_bfai1/2)));
k2_bfai1 = h*(by*(vn+k1_v/2)*zigama1 + (g*sin(theta2n+k1_theta2/2)*sin(bfai1n+k1_bfai1/2)/(vn+k1_v/2)));

k3_x = h*(vn+k2_v/2)*cos(theta2n+k2_theta2/2)*cos(bfai1n+k2_bfai1/2);
k3_y = h*(vn+k2_v/2)*sin(theta2n+k2_theta2/2)*cos(bfai1n+k2_bfai1/2);
k3_z = h*(vn+k2_v/2)*sin(bfai1n+k2_bfai1/2);
k3_v = h*(-bx)*(vn+k2_v/2)^2-g*sin(theta2n+k2_theta2/2)*cos(bfai1n+k2_bfai1/2);
k3_theta2 = h*((by*(vn+k2_v/2)^2*zigama2-g*cos(theta2n+k2_theta2/2))/((vn+k2_v/2)*cos(bfai1n+k2_bfai1/2)));
k3_bfai1 = h*(by*(vn+k2_v/2)*zigama1 + (g*sin(theta2n+k2_theta2/2)*sin(bfai1n+k2_bfai1/2)/(vn+k2_v/2)));

k4_x = h*(vn+k3_v)*cos(theta2n+k3_theta2)*cos(bfai1n);
k4_y = h*(vn+k3_v)*sin(theta2n+k3_theta2)*cos(bfai1n);
k4_z = h*(vn+k3_v)*sin(bfai1n);
k4_v = h*(-bx)*(vn+k3_v)^2-g*sin(theta2n+k3_theta2)*cos(bfai1n);
k4_theta2 = h*((by*(vn+k3_v)^2*zigama2-g*cos(theta2n+k3_theta2))/((vn+k3_v)*cos(bfai1n)));
k4_bfai1 = h*(by*(vn+k3_v)*zigama1 + (g*sin(theta2n+k3_theta2)*sin(bfai1n)/(vn+k3_v)));

xn = xn + (k1_x+2*k2_x+2*k3_x+k4_x)/6;
yn = yn + (k1_y+2*k2_y+2*k3_y+k4_y)/6;
zn = zn + (k1_z+2*k2_z+2*k3_z+k4_z)/6;
vn = vn + (k1_v+2*k2_v+2*k3_v+k4_v)/6;
theta2n = theta2n + (k1_theta2+2*k2_theta2+2*k3_theta2+k4_theta2)/6;
bfai1n = bfai1n + (k1_bfai1+2*k2_bfai1+2*k3_bfai1+k4_bfai1)/6;

xout = [xout,xn];
yout = [yout, yn];
zout = [zout, zn];
vout = [vout, vn];
theta2out = [theta2out, theta2n];
bfai1out = [bfai1out, bfai1n];

end

2.测试程序(因为我不知道每个参数具体代表什么，初值我就随便设的)

clc,clear;
[xout, yout, zout, vout, theta2out, bfai1out] = Runge_Kutta4(10,0.1,0.1,0,0,0,4,0,0,0.1,0.1,0.4,1,10,0.01);
hold on;
plot(xout);
plot(yout);
plot(zout);
plot(vout);
legend('x','y','z','v')