龍格庫塔的算力詳解
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');
具體參數自己設置