展示一下RK6的公式
rk6法误差较大
2个视图（过去30天）
显示旧的评论
下面的代码解决了使用RK6没有阻力的抛射运动。我发现从RK6得到的值和从精确解得到的值之间有很大的差异。有人能解释为什么会这样吗?
清晰的所有;
%用RK6求解无阻力弹丸运动
%x " = 0, y " = -g待解方程
g = 9.807;
韦尔= 100;th_deg = 30;%的输入
x0 = 0;y0 = 0;%初始条件
t0=0；tf=100；%时间跨度
vx =韦尔* cosd (th_deg);沿着x %的速度
v =韦尔*信德(th_deg);%的速度沿着y
将二阶微分方程转化为一阶微分方程
%x'=u用fx表示
%u'=0用fu表示
%y'=v由fy表示
%v'=-g由fv表示
fx=@（t，x，u）u；
傅= @ (t, x, u) 0;
fy=@（t，y，v）v；
阵线= @ (t y v) - g;
t (1) = 0;
x（1）=0；y（1）=0；
x_exact (1) = 0; y_exact (1) = 0;%精确解
u（1）=vx；v（1）=vy；
h = 0.001;
N =装天花板(tf-t (1)) / h);
为j=1:N
t (j + 1) = t (j) + h;
x_exact (j + 1) = vx * t (j + 1);
精确（j+1）=vy*t（j+1）-g*0.5*（t（j+1））^2；
k1x =外汇(t (j) x (j), u (j));
k1u=fu（t（j），x（j），u（j））；
k1y =财政年度(t (j), y (j), v (j));
k1v =阵线(t (j), y (j), v (j));
k2x =外汇(t (j) + h x (j) + k1x u (j) + k1u);
傅k2u = (t (j) + h, x (j) + k1x u (j) + k1u);
k2y =财政年度(t (j) + h, y (j) + k1y v (j) + k1v);
k2v =阵线(t (j) + h, y (j) + k1y v (j) + k1v);
k3x=fx（t（j）+h/2，x（j）+（3*k1x+k2x）/8），u（j）+（3*k1u+k2u）/8）；
k3u=fu（t（j）+h/2，x（j）+（3*k1x+k2x）/8），u（j）+（3*k1u+k2u）/8）；
k3y =财政年度(t (j) + h / 2, y (j) + ((3 * k1y + k2y) / 8), v (j) + ((3 * k1v + k2v) / 8));
k3v =阵线(t (j) + h / 2, y (j) + ((3 * k1y + k2y) / 8), v (j) + ((3 * k1v + k2v) / 8));
k4x =外汇(t (j) + 2 * h / 3、x (j) + (8 * k1x + 2 * k2x + 8 * k3x) / 27日u (j) + (8 * k1u + 2 * k2u + 8 * k3u) / 27);
傅k4u = (t (j) + 2 * h / 3、x (j) + (8 * k1x + 2 * k2x + 8 * k3x) / 27日u (j) + (8 * k1u + 2 * k2u + 8 * k3u) / 27);
k4y =财政年度(t (j) + 2 * h / 3, y (j) + (8 * k1y + 2 * k2y + 8 * k3y) / 27 v (j) + (8 * k1v + 2 * k2v + 8 * k3v) / 27);
k4v=fv（t（j）+2*h/3，y（j）+（8*k1y+2*k2y+8*k3y）/27，v（j）+（8*k1v+2*k2v+8*k3v）/27）；
k5x=fx（t（j）+（7-（21）^0.5）*h/14，x（j）+（3*（3*（21）^0.5-7）*k1x-8*（7-（21）^0.5）*k2x+48*（7-（21）^0.5）*k3x-3*（21-（21）^0.5）*k4x）/392，u（j）+3*（3*（3*（3*（3*（21）^0.5-7*（21）^0.5）/k8u-8*（7-（21）^0.5）*+48*（7-（7-（21）^0.5）/k21）/k2u*0.5*（21）/k21）/k21）；
傅k5u = (t (j) + (7 - (21) ^ 0.5) * h / 14日x (j) + (3 * 3 * 7 (21) ^ 0.5) * k1x-8 * (7 - (21) ^ 0.5) * k2x + 48 * (7 - (21) ^ 0.5) * k3x-3 * (21 - (21) ^ 0.5) * k4x) / 392, u (j) + 3 * ((3 * 7 (21) ^ 0.5) * k1u-8 * (7 - (21) ^ 0.5) * k2u + 48 * (7 - (21) ^ 0.5) * k3u-3 * (21 - (21) ^ 0.5) * k4u) / 392);
k5y=fy（t（j）+（7-（21）^0.5）*h/14，y（j）+（3*（3*（21）^0.5-7）*k1y-8*（7-（21）^0.5）*k2y+48*（7-（21）^0.5）*k3y-3*（21-（21）^0.5）*k4y）/392，v（j）+3*（3*（3*（3*（21）^0.5-7*（7-（21）^0.5）*+48*（7-（7-（21）^0.5）/k3y-5*（21）/k21）/k2*（21）/k21）；
k5v=fv（t（j）+（7-（21）^0.5）*h/14，y（j）+（3*（3*（21）^0.5-7）*k1y-8*（7-（21）^0.5）*k2y+48*（7-（21）^0.5）*k3y-3*（21-（21）^0.5）*k4y）/392，v（j）+3*（3*（3*（3*（21）^0.5-7*（7-（21）^0.5）*+48*（7-（21）^0.5）/k3y*（21）^0.5）*（k21）/k2y）/k21）；
k6x =外汇(t (j) + (7 + (21) ^ 0.5) * h / 14日x (j) + (5 * (231 + 51 * (21) ^ 0.5) * k1x-40 * (7 + (21) ^ 0.5) * k2x - 320 * (21) ^ 0.5 * k3x + 3 * (21 + 121 * (21) ^ 0.5) * k4x + 392 * (6 + (21) ^ 0.5) * k5x) / 1960 u (j) + (5 * (231 + 51 * (21) ^ 0.5) * k1u-40 * (7 + (21) ^ 0.5) * k2u - 320 * (21) ^ 0.5 * k3u + 3 * (21 + 121 * (21) ^ 0.5) * k4u + 392 * (6 + (21) ^ 0.5) * k5u) / 1960);
傅k6u = (t (j) + (7 + (21) ^ 0.5) * h / 14日x (j) + (5 * (231 + 51 * (21) ^ 0.5) * k1x-40 * (7 + (21) ^ 0.5) * k2x - 320 * (21) ^ 0.5 * k3x + 3 * (21 + 121 * (21) ^ 0.5) * k4x + 392 * (6 + (21) ^ 0.5) * k5x) / 1960 u (j) + (5 * (231 + 51 * (21) ^ 0.5) * k1u-40 * (7 + (21) ^ 0.5) * k2u - 320 * (21) ^ 0.5 * k3u + 3 * (21 + 121 * (21) ^ 0.5) * k4u + 392 * (6 + (21) ^ 0.5) * k5u) / 1960);
k6y =财政年度(t (j) + (7 + (21) ^ 0.5) * h / 14日y (j) + (5 * (231 + 51 * (21) ^ 0.5) * k1y-40 * (7 + (21) ^ 0.5) * k2y - 320 * (21) ^ 0.5 * k3y + 3 * (21 + 121 * (21) ^ 0.5) * k4y + 392 * (6 + (21) ^ 0.5) * k5y) / 1960 v (j) + (5 * (231 + 51 * (21) ^ 0.5) * k1v-40 * (7 + (21) ^ 0.5) * k2v - 320 * (21) ^ 0.5 * k3v + 3 * (21 + 121 * (21) ^ 0.5) * k4v + 392 * (6 + (21) ^ 0.5) * k5v) / 1960);
k6v =阵线(t (j) + (7 + (21) ^ 0.5) * h / 14日y (j) + (5 * (231 + 51 * (21) ^ 0.5) * k1y-40 * (7 + (21) ^ 0.5) * k2y - 320 * (21) ^ 0.5 * k3y + 3 * (21 + 121 * (21) ^ 0.5) * k4y + 392 * (6 + (21) ^ 0.5) * k5y) / 1960 v (j) + (5 * (231 + 51 * (21) ^ 0.5) * k1v-40 * (7 + (21) ^ 0.5) * k2v - 320 * (21) ^ 0.5 * k3v + 3 * (21 + 121 * (21) ^ 0.5) * k4v + 392 * (6 + (21) ^ 0.5) * k5v) / 1960);
k7x =外汇(t (j) + h x (j) + (15 * (22 + 7 * (21) ^ 0.5) * k1x + 120 * k2x + 40 * (7 * (21) ^ 0.5 5) * k3x - 63 * (3 * 0.5 (21) ^ 2) * k4x-14 * (49 + 9 * (21) ^ 0.5) * k5x + 70 * (7 - (21) ^ 0.5) * k6x) / 180 u (j) + (15 * (22 + 7 * (21) ^ 0.5) * k1u + 120 * k2u + 40 * (7 * (21) ^ 0.5 5) * k3u - 63 * (3 * 0.5 (21) ^ 2) * k4u-14 * (49 + 9 * (21) ^ 0.5) * k5u + 70 * (7 - (21) ^ 0.5) * k6u) / 180);
傅k7u = (t (j) + h x (j) + (15 * (22 + 7 * (21) ^ 0.5) * k1x + 120 * k2x + 40 * (7 * (21) ^ 0.5 5) * k3x - 63 * (3 * 0.5 (21) ^ 2) * k4x-14 * (49 + 9 * (21) ^ 0.5) * k5x + 70 * (7 - (21) ^ 0.5) * k6x) / 180, u (j) + (15 * (22 + 7 * (21) ^ 0.5) * k1u + 120 * k2u + 40 * (7 * (21) ^ 0.5 5) * k3u - 63 * (3 * 0.5 (21) ^ 2) * k4u-14 * (49 + 9 * (21) ^ 0.5) * k5u + 70 * (7 - (21) ^ 0.5) * k6u) / 180);
k7y =财政年度(t (j) + h, y (j) + (15 * (22 + 7 * (21) ^ 0.5) * k1y + 120 * k2y + 40 * (7 * (21) ^ 0.5 5) * k3y - 63 * (3 * 0.5 (21) ^ 2) * k4y-14 * (49 + 9 * (21) ^ 0.5) * k5y + 70 * (7 - (21) ^ 0.5) * k6y) / 180 v (j) + (15 * (22 + 7 * (21) ^ 0.5) * k1v + 120 * k2v + 40 * (7 * (21) ^ 0.5 5) * k3v - 63 * (3 * 0.5 (21) ^ 2) * k4v-14 * (49 + 9 * (21) ^ 0.5) * k5v + 70 * (7 - (21) ^ 0.5) * k6v) / 180);
k7v =阵线(t (j) + h, y (j) + (15 * (22 + 7 * (21) ^ 0.5) * k1y + 120 * k2y + 40 * (7 * (21) ^ 0.5 5) * k3y - 63 * (3 * 0.5 (21) ^ 2) * k4y-14 * (49 + 9 * (21) ^ 0.5) * k5y + 70 * (7 - (21) ^ 0.5) * k6y) / 180 v (j) + (15 * (22 + 7 * (21) ^ 0.5) * k1v + 120 * k2v + 40 * (7 * (21) ^ 0.5 5) * k3v - 63 * (3 * 0.5 (21) ^ 2) * k4v-14 * (49 + 9 * (21) ^ 0.5) * k5v + 70 * (7 - (21) ^ 0.5) * k6v) / 180);
x（j+1）=x（j）+h*（9*k1x+64*k3x+49*k5x+49*k6x+9*k7x）/180；
u（j+1）=u（j）+h*（9*k1u+64*k3u+49*k5u+49*k6u+9*k7u）/180；
y（j+1）=y（j）+h*（9*k1y+64*k3y+49*k5y+49*k6y+9*k7y）/180；
v（j+1）=v（j）+h*（9*k1v+64*k3v+49*k5v+49*k6v+9*k7v）/180；
如果（y（j+1）<0）
打破;
结束
结束
hmax = max (y);
征求= max (x);
情节(x, y);
持有在;
情节(x_exact y_exact);
包含(“X”)；
ylabel (“Y”)；
答案(2)
darova
2021年8月9日
首先，你应该重写你的代码，使其更具可读性。代码很难检查。尽可能把它讲清楚
t (j + 1) = t (j) + h;
x_exact (j + 1) = vx * t (j + 1);
精确（j+1）=vy*t（j+1）-g*0.5*（t（j+1））^2；
tj = t (j);
xj=x（j）；
yj = y (j);
uj = u (j);
vj=v（j）；
k1x =外汇(tj, xj, uj);
% ... 类比地
k2x =外汇(tj + h, xj + k1x uj + k1u);
%……
C3x = (3 * k1x + k2x) / 8;
C3u=（3*k1u+k2u）/8；
%……另一个C3…
k3x =外汇(tj + h / 2, xj + C3x uj + C3u);
%……
%……
C71 = 15 * (22 + 7 * (21) ^ 0.5);
C73 = 40 * (7 * (21) ^ 0.5 5);
C74 = 63 * 3 * 0.5 (21) ^ 2);
C75 = 14 * (49 + 9 * (21) ^ 0.5);
C76=70*（7-（21）^0.5）；
dx7=（C71*k1x+120*k2x+C73*k3x-C74*k4x-C75*k5x+C76*k6x）/180；
du7 = (C71 * k1u + 120 * k2u + C73 * k3u-C74 * k4u-C75 * k5u + C76 * k6u) / 180;
tj+h, xj+dx7, uj+du7;
%……
2评论
darova
2021年8月11日
我犯不了一个错误。这里是比较
数值
和
确切的
解决方案。看起来算法中有错误
RK6
g = 9.807;
韦尔= 100;th_deg = 30;%的输入
x0 = 0;y0 = 0;%初始条件
t0=0；tf=100；%时间跨度
vx =韦尔* cosd (th_deg);沿着x %的速度
v =韦尔*信德(th_deg);%的速度沿着y
t=linspace（0,9）；
x_exact = vx * t;%精确解x=vx*t
y_精确=vy*t-g*0.5*t.^2；%精确解y=vy*t+0.5*g*t^2
F = @(t, F) [F (3); F (4);0;];
[t1,f1] = ode45(f,t，[0 0 vx vy]);
绘图（x_精确，y_精确，“r”)
线(f1 f1 (: 1), (2):,)
传奇(精确解的,“数值”)
你也可以从以下列表中选择一个网站: