问题标签 [runge-kutta]

有人可以向我解释为什么 Verlet 积分比 Euler 积分更好吗？为什么 RK4 比 Verlet 更好？我不明白为什么这是一个更好的方法。

这个方程主要是这个线程的结果：Java 中的微分方程。
基本上，我尝试遵循 Jason S. 的建议并通过 Runge-Kutta 方法 (RK4) 实现微分方程的数值解。

大家好，我正在尝试用java创建一个简单的SIR-epidemics模型模拟程序。基本上，SIR 由三个微分方程组定义：
S'(t) = - lamda(t) * S(t)
I'(t) = lamda(t) * S(t) - gamma(t) * I(t)
R'(t) = gamma(t) * I(t)
S - 易感人群，I - 感染者，R - 康复者。lamda(t) = [c * x * I(t)] / N(T) c - 接触人数，x - 传染性（与病人接触后生病的概率），N(t) - 总人口（其中是恒定的）。
gamma(t) = 1 / 病程（常数）

在第一次不是很成功的尝试之后，我尝试用 Runge-KUtta 解决这个方程，这个尝试产生了以下代码：

这似乎不起作用，因为生病的人数（I）应该先增加，然后减少到大约0，并且恢复的人数应该严格增加。生病 + 健康 + 康复的总数应该是 100，但我的代码产生了一些奇怪的结果：

找不到错误，请指教！提前谢谢了！

我正在尝试在二维矩阵（x，y）上的 matlab 中实现 Runge-Kutta 的四阶积分方法。该矩阵在每个点都包含一个高度值（浮点数）。这个想法是在矩阵中放置一个粒子并观察它的路径，因为它被方向向量 u 和 v “推动”。

我已经简单地使用 matlab 中内置的“梯度”函数实现了欧拉的积分方法。但是，对于 RK4 积分，我需要在可变点取 4 个斜率，而不是梯度给出的预定斜率。

我认为 Matlab 梯度函数本身没有这些参数（根据我从 Matlab 文档中了解到的）。我应该采取什么样的方法？

如果我要编写自己的梯度函数来计算任意点的斜率，我不太确定如何编写它，因为它不是一个直接的函数。对这种方法有什么建议吗？

谢谢。

我试图确定Runge-Kutta 方法（“RK4”）在常微分方程精确解的 0.01% 范围内需要多少步。我将此与欧拉方法进行比较。两者都应该在对数图上产生一条直线。我的欧拉解决方案似乎是正确的，但我得到了 RK 的弯曲解决方案。它们基于相同的代码，所以我对这个问题完全感到困惑。

编辑：很抱歉删除了维基百科链接。由于我是新用户，它不会让我保留多个链接。但是，这两种方法在维基百科上都有详细介绍，类似于我的实现。

如果有人希望解决我的问题，代码如下，图表位于dropbox.com 上的 Word 文件中。是的，这是一个家庭作业问题；我发布这个是因为我想了解我的思考过程中出了什么问题。

这是我的 RK4 代码：

这是我的欧拉代码：

我正在尝试使用 Runge Kutta 方法计算 Lorenz 系统，但我找不到我的代码哪里有错误。当我运行它时，系统会进入一个静态点，我应该获得一只蝴蝶（洛伦兹吸引子）。我认为这是 Runge Kutta 方法部分的“for”循环中的内容。这是我的代码

提前致谢

编辑

我在 linux 上再次（以更简化的方式）编写了代码，它运行正常，似乎我在 Windows 上的编辑器有一些东西（可能是编码），当我编译代码时，它抛出了很多无穷大的值。不知道为什么，我花了很多时间才注意到这一点。谢谢你的帮助

这是我的第一篇文章，我不得不承认，我在编程方面很糟糕。我是班上那个拼命工作的人，但似乎永远无法像其他同学一样掌握编程。所以请友善，我将在下面尝试解释我的问题。

我有以下代码（已删除注释），但是当我运行它时，我收到类似于下面列出的警告。另外，当我运行程序时，允许第一个用户输入的值，但突然之间，它跳到程序的末尾，不允许我输入其他变量的值（例如变量“beta” ）。我有输出图像（http://i.stack.imgur.com/yc3jq.jpg），您可以看到我输入了 alpha，但随后程序运行到最后。有什么想法吗？

非常感谢您的帮助！-斯宾塞

- - - - - - - - - - - - - - -代码 - - - - - - - -

示例警告：

有谁知道使用 runge kutta 的 2D 刚体物理的任何 c++/opengl 源代码演示？

我想构建一个物理引擎，但我需要一些参考代码来更好地理解其他人是如何实现的。

自从我做任何 C/c++ 以来已经有一段时间了，但我想使用 gsl 库编写一个 ODE 求解器来求解以下 ODE 集

所以在 gsl 表示法中我的 y[0]=u, y[1]==up，上面的 RHS 定义了 f[0] 和 f[1]。根据这些定义，可以计算雅可比和 dfdr（通常它们的“时间”变量称为“t”而不是“r”）。这样做的原因是因为我在使用 Mathematica 时遇到了速度问题。我在他们关于 ODE 求解器的文档末尾获取了 gsl 示例代码，并尝试使其适应我的问题，如下所示：

这是给出的数字，但它们与 Mathematica NDSolve 给出的数字不同，即使我的 WorkingPrecision 和 PrecisionGoal 较低。我所做的是否有错误？

根据在这个答案中给我的建议，我在我的重力模拟器中实现了一个龙格-库塔积分器。

然而，在我模拟了太阳系一年之后，位置仍然偏离了 110 000 公里，这是不可接受的。

我的初始数据由 NASA 的 HORIZONS 系统提供。通过它，我获得了行星、冥王星、月球、火卫二和火卫一在特定时间点的位置和速度矢量。

这些矢量是 3D 的，但是，有些人告诉我，我可以忽略三维，因为行星在围绕太阳的盘子中排列自己，所以我做到了。我只是将 xy 坐标复制到我的文件中。

这是我改进的更新方法的代码：

算法如下：

1. 更新系统中所有物体的加速度
2. RK4(第一步)
3. 转到 1
4. RK4（二）
5. 转到 1
6. RK4（第三）
7. 转到 1
8. RK4（第四）

我的 RK4 实施是否搞砸了？还是我只是从损坏的数据开始（重要的机构太少而忽略了第三维）？

如何解决这个问题？

解释我的数据等...

我所有的坐标都是相对于太阳的（即太阳在 (0, 0)）。

我110 000 km通过从我的模拟器预测的坐标中减去 NASA 给出的地球 x 坐标得到了错误。

相对误差似乎微乎其微，但这仅仅是因为在我的模拟和美国宇航局的模拟中地球都离太阳很远。距离仍然很大，使我的模拟器毫无用处。

假设我有一条贝塞尔曲线 B(u)，如果我u以恒定速率增加参数，我不会沿着曲线获得恒定的速度运动，因为u参数和获得的评估曲线的点之间的关系不是线性的。

我已经阅读并实施了 David Eberly 的文章。它解释了如何沿着参数曲线以恒定速度移动。

假设我有一个F(t)将时间值作为输入的函数和一个返回时间值t的速度函数，我可以获得沿曲线的恒定速度运动，以恒定速率改变 t 参数：sigmatB(F(t))

我正在使用的文章的核心是以下功能：

它允许我u使用提供的时间t和 sigma 函数计算曲线参数。现在，当速度 sigma 为常数时，该函数可以正常工作。如果 sigma 代表一个统一的加速，我会从中得到错误的值。

这是一个直线贝塞尔曲线的示例，其中 P0 和 P1 是控制点，T0 T1 是切线。曲线的定义：

假设我想知道 time 沿曲线的位置t = 3。如果我是一个恒定的速度：

以及以下数据：

我可以分析计算位置：

如果我使用我的贝塞尔样条曲线和上面的算法求解相同的方程，n =5我得到：

考虑到数值近似，该值非常精确（我对此做了很多测试。我省略了细节，但贝塞尔曲线实现很好，曲线本身的长度是使用Gaussian Quadrature非常精确地计算出来的）。

现在，如果我尝试将 sigma 定义为均匀加速函数，我会得到不好的结果。考虑以下数据：

我可以使用线性运动方程计算粒子到达 P1 的时间：

t我可以计算加速度：

我有所有数据来定义我的 sigma 函数：

如果我分析地解决这个问题，我希望在 time 之后粒子的以下速度t =3：

职位将是：

但如果我用上面的算法计算它，位置结果：

这与我的预期完全不同。

很抱歉这篇文章很长，如果有人有足够的耐心阅读它，我会很高兴。

你有什么建议吗？我错过了什么？谁能告诉我我做错了什么？

编辑：我正在尝试使用 3D 贝塞尔曲线。这样定义：

和导数：