Cómo utilizar el método Runge-Kutta para resolver problemas mecánicos, gracias. La pregunta está en la imagen. Se requieren procedimientos detallados, incluido el dibujo del diagrama.
Para utilizar el método RK, primero debe obtener la ecuación diferencial ordinaria del problema. Dado que el problema es muy simple, el análisis y la solución específicos no se discutirán aquí, es relevante para esto. problema
v'= g-k*v^2/m
Entre ellos, v' es la derivada de la velocidad con respecto al tiempo, es decir, la aceleración
Usando el método RK4, dado que la ecuación diferencial ordinaria solo está relacionada con v, está relacionada con t. No existe una relación directa. El método RK4 se puede escribir simplemente como
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)
Según la pregunta, hay un valor inicial v(0)=0 (la velocidad inicial es 0), eso es todo Solución
Como no tengo VB instalado a mano, Reescribí el programa en lenguaje Tcl. Tcl / tk puede ejecutarse normalmente. De todos modos, puede consultarlo. Las ideas de programación son similares.
Adjunte primero el escrito en Tcl/Tk:
set g 9.81
set k 29.4
set m 75
#Definir funciones diferenciales ordinarias
Proc f {v} {
global g k m
set f [expr $g-$k*$v * *2/$m]
}
#La siguiente es la parte principal del cálculo, llamando a la función anterior
set h 0.1
establecer v(0) 0
establecer n 0
establecer dv 1
mientras {$dv > 10**-3} {
establecer k1 [f $v($n)]
establecer k2 [f [expr $v($n) + $h / 2 * $k1]]
establecer k3 [f [expr $v($n) + $h / 2 * $k2]]
establecer k4 [f [expr $v($n) + $h * $k3]]
establecer nn [expr $n+1]
establecer v($nn) [expr $v($n) + $h / 6 * ($k1 + 2 * $ k2 + 2 * $k3 + $k4)]
pone [expr $h * $n],$v($nn)
establece dv [expr abs($v( $n ) - $v($nn))]
incr n
}
Resultado de ejecución de Tcl/Tk:
Tiempo , Velocidad
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 32
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 85401
2.8 5.0023797583
2.9 5.002435097
p>
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,0 882609
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,002550 34929
5.2 5.002550356
5.3 5.00255036053
5.4 5.00255036359
5.5 5.00255036566
5.6 5.00255036706
5,7 5,0 02550368
5,8 5,00255036864
5,9 5,00255036907
6,0 5,00255036936
6,1 5,00255036956
6,2 36969
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
La solución exacta es la velocidad máxima 5.002550369969
El siguiente es el código VB (no depurado):
Const privada g = 9,81
Const privada k = 29,4
Const privada m = 75
Función f(ByVal v) siempre
f = g - k * v ^ 2 / m
Función final
Sub privada RK4()
Atenuar h, k1, k2, k3, k4, dv siempre
Atenuar v(100) mientras
h = 0,1
v( 0) = 0
n = 0
dv = 1
Hacer mientras d
v > 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
Bucle
End Sub
Creo que es sencillo Para producir la imagen es más fácil escribir en Matlab o dibujar con Excel, así que no escribiré el código relevante. Después de todo, esa parte es mucho más complicada.