龙格库塔的算力详解
A. 龙格库塔法的计算公式
（1）的局部截断误差是 。
龙格-库塔法具有精度高，收敛，稳定（在一定条件下），计算过程中可以改变步长，不需要计算高阶导数等优点，但仍需计算 在一些点上的值，如四阶龙格-库塔法每计算一步需要计算四次 的值，这给实际计算带来一定的复杂性，因此，多用来计算“表头”。
二、小程序
#include<stdio.h>
#include<math.h>
#define f(x,y) (-1*(x)*(y)*(y))
void main(void)
{
double a,b,x0,y0,k1,k2,k3,k4,h;
int n,i;
printf(input a,b,x0,y0,n:);
scanf(%lf%lf%lf%lf%d,&a,&b,&x0,&y0,&n);
printf(x0 y0 k1 k2 k3 k4
);
for(h=(b-a)/n,i=0;i!=n;i++)
{
k1=f(x0,y0);
k2=f(x0+h/2,y0+k1*h/2);
k3=f(x0+h/2,y0+k2*h/2);
k4=f(x0+h,y0+h*k3);
printf(%lf %lf ,x0,y0);
printf(%lf %lf ,k1,k2);
printf(%lf %lf
,k3,k4);
y0+=h*(k1+2*k2+2*k3+k4)/6;
x0+=h;
}
printf(xn=%lf yn=%lf
,x0,y0);
}
运行结果：
input a,b,x0,y0,n:0 5 0 2 20
x0 y0 k1 k2 k3 k4
0.000000 2.000000 -0.000000 -0.500000 -0.469238
-0.886131
0.250000 1.882308 -0.885771 -1.176945 -1.129082
-1.280060
0.500000 1.599896 -1.279834 -1.295851 -1.292250
-1.222728
0.750000 1.279948 -1.228700 -1.110102 -1.139515
-0.990162
1.000000 1.000027 -1.000054 -0.861368 -0.895837
-0.752852
1.250000 0.780556 -0.761584 -0.645858 -0.673410
-0.562189
1.500000 0.615459 -0.568185 -0.481668 -0.500993
-0.420537
1.750000 0.492374 -0.424257 -0.361915 -0.374868
-0.317855
2.000000 0.400054 -0.320087 -0.275466 -0.284067
-0.243598
2.250000 0.329940 -0.244935 -0.212786 -0.218538
-0.189482
2.500000 0.275895 -0.190295 -0.166841 -0.170744
-0.149563
2.750000 0.233602 -0.150068 -0.132704 -0.135399
-0.119703
3.000000 0.200020 -0.120024 -0.106973 -0.108868
-0.097048
3.250000 0.172989 -0.097256 -0.087300 -0.088657
-0.079618
3.500000 0.150956 -0.079757 -0.072054 -0.073042
-0.066030
3.750000 0.132790 -0.066124 -0.060087 -0.060818
-0.055305
4.000000 0.117655 -0.055371 -0.050580 -0.051129
-0.046743
4.250000 0.104924 -0.046789 -0.042945 -0.043363
-0.039833
4.500000 0.094123 -0.039866 -0.036750 -0.037072
-0.034202
4.750000 0.084885 -0.034226 -0.031675 -0.031926
-0.029571
xn=5.000000 yn=0.076927
B. 龙格库塔求解微分方程组
只要理解了龙格库塔，这就很容易了。
定义函数，
f=func(x,y)
if y(1)>126,
f=[f1(y(1),y(2)), f2(y(1),y(2))];
else
f=[g1(y(1),y(2)), g2(y(1),y(2))];
end
然后就简单了
C. 跪求四阶龙格库塔公式求 常微分方程的方法 不要代码 只要具体的推倒或者是 解出的方程组
http://ke..com/link?url=dTC8pz6PJXWpZdsua4rFN_panquCZBoHgwqPNwdvVqgZL_aAEXaCy7vebNX-Sk_-rPmjGGne0mNuaYTurUhqP90-_qfKXnpwaxfKcyUt
就有推导啊
指出楼主的一个认识错误 龙格库塔法仅仅是数值求解 并不能说是解出了方程组 数值解再怎么说也只是解析解的近似描述 （个人认为 仅供参考）
D. 龙格库塔算法是什么
参考答案 别在喜悦时许诺，别在忧伤时回答，别在愤怒时做决定。
E. 怎么用龙格库塔法求解力学问题，谢谢，题目在图片里，要详细的程序，包括画出图
要使用RK法，首先要得到问题的常微分方程，具体的分析求解因为题目很简单，这里就不说了，反正本题有关系
v'=g-k*v^2/m
其中，v'为速度对时间的导数，也就是加速度
采用RK4法，由于常微分方程只与v有关，与t无直接关系，RK4法可以简单写为
k1=f(v(n))
k2=f(v(n)+h/2*k1)
k3=f(v(n)+h/2*k2)
k4=f(v(n)+h*k3)
v(n+1)=v(n)+h/6*(k1+2k2+2k3+k4)
根据题目，有初值v(0)=0 (初速度为0），即可求解
由于手头没有装VB，所以我把一下程序用Tcl语言重写了一遍，Tcl/tk的可以正常运行，反正你可以参照着看，所有的语言编程思路是相似的。
先附上Tcl/Tk写的：
set g 9.81
set k 29.4
set m 75
#定义常微分函数
Proc f {v} {
global g k m
set f [expr $g-$k*$v**2/$m]
}
#以下是主要计算部分，调用了上面的函数
set h 0.1
set v(0) 0
set n 0
set dv 1
while {$dv > 10**-3} {
set k1 [f $v($n)]
set k2 [f [expr $v($n) + $h / 2 * $k1]]
set k3 [f [expr $v($n) + $h / 2 * $k2]]
set k4 [f [expr $v($n) + $h * $k3]]
set nn [expr $n+1]
set v($nn) [expr $v($n) + $h / 6 * ($k1 + 2 * $k2 + 2 * $k3 + $k4)]
puts [expr $h * $n],$v($nn)
set dv [expr abs($v($n) - $v($nn))]
incr n
}
Tcl/Tk的运行结果：
时间，速度
0.1 0.968603324806
0.2 1.86719840948
0.3 2.64465167179
0.4 3.27770943734
0.5 3.76821957195
0.6 4.13386114838
0.7 4.39864064439
0.8 4.58638649863
0.9 4.71753254441
1.0 4.80818704498
1.1 4.87039891085
1.2 4.9128797342
1.3 4.94178878159
1.4 4.96141640673
1.5 4.97472149869
1.6 4.98373107432
1.7 4.98982752078
1.8 4.99395074243
1.9 4.99673848624
2.0 4.99862288083
2.1 4.99989645764
2.2 5.00075712243
2.3 5.00133870705
2.4 5.00173168802
2.5 5.00199721975
2.6 5.002176632
2.7 5.00229785401
2.8 5.0023797583
2.9 5.002435097
3.0 5.00247248647
3.1 5.00249774851
3.2 5.00251481666
3.3 5.00252634865
3.4 5.00253414015
3.5 5.00253940442
3.6 5.00254296119
3.7 5.00254536429
3.8 5.00254698792
3.9 5.00254808492
4.0 5.00254882609
4.1 5.00254932686
4.2 5.0025496652
4.3 5.0025498938
4.4 5.00255004825
4.5 5.0025501526
4.6 5.0025502231
4.7 5.00255027074
4.8 5.00255030293
4.9 5.00255032467
5.0 5.00255033936
5.1 5.00255034929
5.2 5.002550356
5.3 5.00255036053
5.4 5.00255036359
5.5 5.00255036566
5.6 5.00255036706
5.7 5.002550368
5.8 5.00255036864
5.9 5.00255036907
6.0 5.00255036936
6.1 5.00255036956
6.2 5.00255036969
6.3 5.00255036978
6.4 5.00255036984
6.5 5.00255036988
6.6 5.00255036991
6.7 5.00255036993
6.8 5.00255036994
6.9 5.00255036995
7.0 5.00255036996
7.1 5.00255036996
精确解为最大速度5.002550369969
以下是VB的代码（未调试）：
Private Const g = 9.81
Private Const k = 29.4
Private Const m = 75
Function f(ByVal v) As Long
f = g - k * v ^ 2 / m
End Function
Private Sub RK4()
Dim h, k1, k2, k3, k4, dv As Long
Dim v(100) As Long
h = 0.1
v(0) = 0
n = 0
dv = 1
Do While dv > 10 ^ -5
k1 = f(v(n))
k2 = f(v(n) + h / 2 * k1)
k3 = f(v(n) + h / 2 * k2)
k4 = f(v(n) + h * k3)
v(n + 1) = v(n) + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
print v(n + 1)
dv = Abs(v(n) - v(n + 1))
n = n + 1
Loop
End Sub
出图我觉得直接在Matlab中编写或者用Excel绘制较简单，所以就不写了相关代码了，毕竟那部分复杂多了
F. 关于二阶和四阶龙格库塔法的计算与编程调试
求解二阶微分方程，初始条件还需要给出y1'(0)和y2'(0)。这里暂时按照0处理。
function zd530003514
a=0.1;
b=0.1;
Y0 = [b-1; 0; b; 0];
% 解方程
[t,Y]= ode45(@ode,[0 10],Y0);
y1=Y(:,1);
y2=Y(:,3);
% 绘图
subplot 211
plot(t,y1);
subplot 212
plot(t,y2);
% 微分方程定义
function dY = ode(t, Y)
L1=5;
L2=0.01;
a0=2;
b0=2;
c0=2;
y1=Y(1);y2=Y(3);
dY = [
Y(2);
-(a0*y2+b0*y2^2+c0*y2^3) - L1^2*L2*y1 - L1^2*y1;
Y(4);
-(a0*y2+b0*y2^2+c0*y2^3) - L1^2*L2*y1;
];
G. MATLAB 龙格库塔法求解常微分方程初值问题 在线等，用MATLAB 写出步骤就行啦~~~！！！！！！ 急
该问题可以求解，可得解析解。dy/dx=a(y-x)+1d(y-x)/dx=a(y-x)即d(y-x)/(y-x)=adx积分得：y-x=Ce^{ax}即：y=Ce^{ax}+x带入初值条件，有C=1从而得特解为：y=e^{ax}+x
H. 龙格库塔三阶公式如何推导
在区间[xi, xi+1]上除两点xi和xi+a2= xi +a2h,以外，再增加一点xi+a3= xi +a3h ，用这三点处的斜率值K1、K2和K3的加权平均得出平均斜率K*的近似值K，这时计算格式具有形式:
y(i+1)=y(i)+c1K1+c2K2+c3K3
K1=hf(xi,yi)
K2=hf(xi+a2h,yi+b21K1)
K3=hf(xi+a3h,yi+b31K1+b32K2)
同理推导二阶公式，将y(xi+1)和yi+1在x=xi处进行Taylor展开，使局部截断误差达到O(h4)，使对应项的系数相等...
I. 用四阶龙格库塔法求解
假设求解初值问题：
y'=y-2x/y
(0<x<1)
y(0)=1
设步长h=0.2,从x=0直到x=1用四阶龙格库塔法：
private
sub
form_click()
dim
x
as
single,
y
as
single
dim
k1
as
single,
k2
as
single,
k3
as
single,
k4
as
single
dim
y1
as
single,
y2
as
single,
y3
as
single,
y4
as
single
x
=
0
y
=
1
h
=
0.2
for
i
=
0
to
1
step
0.2
k1
=
f(x,
y)
y1
=
y
+
h
/
2
*
k1
x
=
x
+
h
/
2
k2
=
f(x,
y1)
y2
=
y
+
h
/
2
*
k2
k3
=
f(x,
y2)
y3
=
y
+
h
*
k3
x
=
x
+
h
/
2
k4
=
f(x,
y3)
y
=
y
+
h
/
6
*
(k1
+
2
*
k2
+
2
*
k3
+
k4)
print
"x=";
x;
"
y=";
y
next
end
sub
private
function
f(a
as
single,
b
as
single)
as
single
f
=
b
-
2
*
a
/
b
end
function
J. 龙格库塔法求解微分方程，matlab怎么编程
function [Y] = RK45(t,X,f,h)
K1=f(t,X);
K2=f(t+h/2,X+h/2*K1);
K3=f(t+h/2,X+h/2*K2);
K4=f(t+h,X+h*K3);
Y=X+h/6*(K1+2*K2+2*K3+K4);
end
以上是4阶龙格库塔法的代码:
自己写函数，存为f.m
function dxdt = f (t,x)
dxdt(1)=exp(x(1)*sin(t))+x(2);
dxdt(2)=exp(x(2)*cos(t))+x(1); % x(1)是你的f，x(2)是你的g
dxdt=dxdt(:);
end
自己给出t0,x0,h的值（初始时间，初值，步长）
如果求t0到t1的轨迹的话:给个例子如下
t0=0;t1=5;h=0.02;x0=[-1;-1];
T=t0:h:t1;X=zeros(length(x0),length(T));X(:,1)=x0;
for j=1:length(T)-1
X(:,j+1)=RK45(T(j),X(:,j),@(t,x) f(t,x),h);
end
plot(T,X(1,:));
hold on;
plot(T,X(2,:),'r');
具体参数自己设置