Question

急急急！求助matlab用龙格-库塔方法求解方程组

微分方程组在图片中，请详细点，谢谢~
麻烦帮忙定义 一下微分方程组，我加分

Answer 1

定义你的微分方程到一个m文件中，然后直接ODE45就可以了
如果是课程作业的话，可能要求你自己写runge-kutta的程序
另外，不知道你所说的RK是几阶，列出几个供你参考（不过不是针对向量运算的，需要自己修改一下）

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% 常微分方程数值解法 (Runge-Kutta)%%%%%%%%%%%%%%%%%%%%%
%%%%% 微分方程Dy=y-2*x/y; y(0)=1 %%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 程序开始，定义求解区间与步长
clear;clc;
h=0.1;% 步长
a=0;% 求解区间下限
b=1;% 求解区间上限
x=a:h:b;% 变量取值
y_analy=sqrt(1+2*x);% 解析解
plot(x,y_analy,'b');hold on; %绘制精确解曲线
n=(b-a)/h;
y0=1;

%% Runge-Kutta 2
% 二阶龙格-库塔法
y_rk2(1)=y0;
for i=1:n
rk2_k1=myode(x(i),y_rk2(i));
rk2_k2=myode(x(i)+h/2,y_rk2(i)+(h/2)*rk2_k1);
y_rk2(i+1)=y_rk2(i)+h*rk2_k2;
end
plot(x,y_rk2,'r')

%% Runge-Kutta 3
% 三阶龙格-库塔法
y_rk3(1)=y0;
for i=1:n
rk3_k1=myode(x(i),y_rk3(i));
rk3_k2=myode(x(i)+h/2,y_rk3(i)+(h/2)*rk3_k1);
rk3_k3=myode(x(i)+h,y_rk3(i)-h*rk3_k1+2*h*rk3_k2);
y_rk3(i+1)=y_rk3(i)+(h/6)*(rk3_k1+4*rk3_k2+rk3_k3);
end
plot(x,y_rk3,'g')

%% Runge-Kutta 4
% 四阶龙格-库塔法
y_rk4(1)=y0;
for i=1:n
rk4_K1=myode(x(i),y_rk4(i));
rk4_K2=myode(x(i)+h/2,y_rk4(i)+(h/2)*rk4_K1);
rk4_K3=myode(x(i)+h/2,y_rk4(i)+(h/2)*rk4_K2);
rk4_K4=myode(x(i)+h,y_rk4(i)+h*rk4_K3);
y_rk4(i+1)=y_rk4(i)+(h/6)*(rk4_K1+2*rk4_K2+2*rk4_K3+rk4_K4);
end
plot(x,y_rk4,'k')

%% Runge-Kutta c
% 变步长龙格-库塔法

% plot(x,y_euler3,'y')

%% analyze error
% 误差分析
err_rk2=y_rk2-y_analy;
err_rk3=y_rk3-y_analy;
err_rk4=y_rk4-y_analy;
% err_ec=y_rkc-y_analy;

figure;
subplot(3,1,1);plot(x,err_rk2)
subplot(3,1,2);plot(x,err_rk3)
subplot(3,1,3);plot(x,err_rk4)
% subplot(4,1,4);plot(x,err_ec)

Answer 2

function df=ode45_fun(t,xyzuvw)
%%注意小写的v和大写的V
%常数(请修正)
R_0=1;
rho_0=1;
beta=1;
G=6.67e-11;
M=1.5e24;

x=xyzuvw(1);
y=xyzuvw(2);
z=xyzuvw(3);
u=xyzuvw(4);
v=xyzuvw(5);
w=xyzuvw(6);

R=sqrt(x*x+y*y+z*z);
V=sqrt(u*u+v*v+w*w);
rho=rho_0*exp(-a*(R-R_0));

first=-rho*V/2/beta;
second=G*M/R^3;

df=[
u;
v;
w;
first*u-second*x;
first*v-second*y;
first*w-second*z;
];