更新了Runge-Kutta(RK4)C++错误代码中的二阶DE

Question

这是Runge Kutta(RK4(在C++错误代码中的二阶DE的更新版本

我仍然在代码方面遇到困难。也许这与我对龙格库塔的有限了解有关，但当我运行这段代码时，它不会产生输出。

#include <iostream> 
#include <cmath>
//dvdt=-(g/L)*sin(theta)
//v=dxdt
double dxdt( double timepassed, double theta )
{
return theta/timepassed;
}
double L;   
double g=9.8;
double coeff=-1*(g/L);
double dvdt( double timepassed,  double x, double v)
{
return coeff*sin(x);
}
int main(){
// Standard Variables   
double theta;
double theta1;
double h = 0.1;
double L;
double timepassed;
double time1;
// Input and Output Printing
std::cout << "Please input initial angle (in decimal radians), length of the pendulum (in meters) and the time desired (in seconds). Click ENTER key after each value."<<"n";
std::cin >> theta1;
std::cin >> L;  
std::cin >> timepassed;
// Specific Variable Declarations
double coeff=-1*(g/L);
double v = dxdt(theta1, timepassed);
double x = theta1;
double d2xdt2 = dvdt(timepassed, theta1, v);
// Defining K Values in Runge Kutta
double kx1,kv1;
double kx2, kv2;
double kx3, kv3;
double kx4, kv4;
double dt;
kx1=dt*dxdt(timepassed,x);
kv1=dt*dvdt(timepassed,x,v);
kx2=dt*dxdt(timepassed+dt/2,x+kx1/2);
kv2=dt*dvdt(timepassed+dt/2,x+kx1/2,v+kv1/2);
kx3=dt*dxdt(timepassed+dt/2,x+kx2/2);
kv3=dt*dvdt(timepassed+dt/2,x+kx2/2,v+kv2/2);
kx4=dt*dxdt(timepassed+dt,x+kx3);
kv4=dt*dvdt(timepassed+dt,x+kx3,v+kv3);
x = x + (1.0/6.0)*(kx1 + 2*kx2 + 2*kx3 + kx4);
v = v + (1.0/6.0)*(kx1 + 2*kv2 + 2*kv3 + kv4);

std::cout << "The angle is" << x; "n"; 
std::cout << "The velocity is" << v;

}

Answer 1

正如前面注释中所宣布的，您的系统方程应该是

//v=dx/dt
//dv/dt=d2x/dt2=-(g/L)*sin(x), where x=theta
double coeff;
double dxdt( double t, double x, double v) { return v; }
double dvdt( double t, double x, double v) { return coeff*sin(x); }

在输入参数之后，数值coeff被计算，但不被重新声明。

// Specific Variable Declarations
coeff = -(g/L);

您的步长显示为0.1。您需要决定使用哪个变量名，h或dt，然后使用它。

几乎可以肯定的是，您需要执行多个RK4步骤，因此您需要用循环来构建它们。除了RK4阶段之外，该循环还包含什么取决于程序的输出。此外，如果目标时间不是时间步长的倍数，则需要调整最后一步。

while(t < timepassed) {
kx1=dt*dxdt(t,x,v);
kv1=dt*dvdt(t,x,v);
kx2=dt*dxdt(t+dt/2,x+kx1/2,v+kv1/2);
kv2=dt*dvdt(t+dt/2,x+kx1/2,v+kv1/2);
kx3=dt*dxdt(t+dt/2,x+kx2/2,v+kv2/2);
kv3=dt*dvdt(t+dt/2,x+kx2/2,v+kv2/2);
kx4=dt*dxdt(t+dt,x+kx3,v+kv3);
kv4=dt*dvdt(t+dt,x+kx3,v+kv3);
t = t + dt;
x = x + (1.0/6.0)*(kx1 + 2*kx2 + 2*kx3 + kx4);
v = v + (1.0/6.0)*(kx1 + 2*kv2 + 2*kv3 + kv4);
// output for the step? Depending on the time?
}