龙格-库塔四阶方法迭代时的累积误差
Runge-Kutta 4th order method cumulative errors when iterating
我试图得到一个简单追逐问题的数值解
(移动目标+带有等速模块的火箭)
每次迭代,我的速度模块都会减小一点,将误差加起来;在几百次迭代之后,错误会激增,速度会急剧下降。
然而,这不是欧拉方法的情况(代码下面的大块),它只在使用RK4方法时弹出。
我不确定错误在哪里以及为什么会发生,所以任何输入都是赞赏的
#include <fstream>
#include <iostream>
#include <vector>
#include <cmath>
#include <cstring>
#define vx(t,x,y) n*V*((t)*(V)-(x))/pow(((t)*(V)-(x))*((t)*(V)-(x))+((h)-(y))*((h)-(y)),0.5)
#define vy(t,y,x) n*V*((h)-(y))/pow(((t)*(V)-(x))*((t)*(V)-(x))+((h)-(y))*((h)-(y)),0.5)
using namespace std;
class Vector {
public:
double x,y;
Vector(double xx, double yy):x(xx), y(yy){};
virtual ~Vector(){}
Vector operator-() {return Vector(-x,-y);};
friend Vector operator-(const Vector &, const Vector &);
friend Vector operator+(const Vector &, const Vector &);
Vector operator*(double l){return Vector(x*l,y*l);};
friend Vector operator*(double, const Vector &);
Vector operator/(double l){return Vector(x/l,y/l);};
void operator+=(const Vector & v ){ x+=v.x; y+=v.y;};
void operator-=(const Vector & v ){ x-=v.x; y-=v.y;};
void operator/=(const Vector & v ){ x/=v.x; y/=v.y;};
friend ostream & operator<<(ostream & os,const Vector & v){os<<"("<<v.x<<", "<<v.y<<")";return os;};
double norm() {return sqrt(x*x+y*y);};
};
Vector operator-(const Vector & v1, const Vector & v2){
return Vector(v1.x-v2.x,v1.y-v2.y);
}
Vector operator+(const Vector & v1, const Vector & v2){
return Vector(v1.x+v2.x,v1.y+v2.y);
}
Vector operator*(double l, const Vector & v){
return Vector(v.x*l,v.y*l);
}
int main() {
Vector posP(0,0);
double V=100.,t = 0,dt = pow(10.,-2),vx,vy,h=1000.,x,y,n=2.,v;
double kx1,kx2,kx3,kx4,ky1,ky2,ky3,ky4;
Vector posE(0,h);
FILE *fp;
fp = fopen("/Users/Philipp/Desktop/a.dat","w");
while(posP.y<(h)){
posE.x=posE.x+V*dt;
x=posP.x;y=posP.y;
kx1 = vx(t,x,y);
ky1 = vy(t,y,x);
kx2 = vx(t+dt/2.0,x+kx1/2.0,y+ky1/2.0);
ky2 = vy(t+dt/2.0,y+ky1/2.0,x+kx1/2.0);
kx3 = vx(t+dt/2.0,x+kx2/2.0,y+ky2/2.0);
ky3 = vy(t+dt/2.0,y+ky2/2.0,x+kx2/2.0);
kx4 = vx(t+dt,x+kx3,y+ky3);
ky4 = vy(t+dt,y+ky3,x+kx3);
posP.x = posP.x + dt*((kx1 + 2.0*(kx2+kx3) + kx4)/6.0);
posP.y = posP.y + dt*((ky1 + 2.0*(ky2+ky3) + ky4)/6.0);
v=sqrt(((kx1 + 2.0*(kx2+kx3) + kx4)/(6.0))*((kx1 + 2.0*(kx2+kx3) + kx4)/(6.0))+((ky1 + 2.0*(ky2+ky3) + ky4)/(6.0))*((ky1 + 2.0*(ky2+ky3) + ky4)/(6.0)));
t+=dt;
if ((posE-posP).norm()<1) break;
fprintf(fp,"%lf %lf %lf %lf n",posP.x, posP.y, v, t);
}
fclose(fp);
return 0;
}
欧拉方法//Euler cycle
while(posP.y<(h)) {
posE.x=posE.x+V*dt;
x=posP.x;y=posP.y;
vx=vx(t,x,y);
vy=vy(t,y,x);
posP.x=posP.x+vx*dt;
posP.y=posP.y+vy*dt;
t+=dt;
if ((posE-posP).norm()<0.1) break;
fprintf(fp,"%lf %lf %lf n",posP.x, posP.y,vx*vx+vy*vy, t);
//速度模块在第三列,你可以看到所有的东西都是200,而不是RK4的情况,甚至第一次迭代它就下降到~199.99985
}
使用
kx2 = vx(t+dt/2.0,x+kx1/2.0,y+ky1/2.0);
但应该是:
kx2 = vx(t+dt/2.0,x+kx1/2.0*dt,y+ky1/2.0*dt);
和后面类似的。或者你可以将所有k值乘以dt:
kx2 = dt*vx(t+dt/2.0,x+kx1/2.0,y+ky1/2.0);
这两种变体对于隐式龙格-库塔方法更加重要
相关文章:
- 使用std::multimap迭代器创建std::list
- 来自 std::list 的迭代器 .end() 按预期返回"0xcdcdcdcdcdcdcdcd"但 .begin()
- C++中带有List类的迭代器Segfault
- 迭代时从向量和内存中删除对象
- 如何在c++迭代器类型中包装std::chrono
- 带过滤器的现代迭代c++集合
- 在c++中检查长方体是否尽可能快地重叠(无迭代)
- C++矢量迭代
- 集合上的输出迭代器:assign和increment迭代器
- Boost Spirit,获取迭代器内部语义动作
- 擦除while循环中迭代的元素
- 实现一个在集合上迭代的模板函数
- 对于set上的循环-获取next元素迭代器
- 在向量内的向量上迭代
- 为什么output_editor Concept不需要output_e迭代器标记
- TSP递归解的迭代形式
- n 维 4 阶龙格-库塔求解器大迭代中的误差
- 在将迭代器与C 中的null进行比较时,编译误差
- 臂组件:循环只会迭代一次,然后较低误差
- 龙格-库塔四阶方法迭代时的累积误差