龙格-库塔四阶
Runge-Kutta 4th-Order
现在我正在努力让我的 Runge-Kutta 发挥作用。实际上我必须解决一个比这个 ODE 系统复杂得多的 ODE 系统,但是 runge-kutta-function 是一样的。我正在解决这个简单的系统,以检查模拟是否收敛到解析解。
首先,如果您发现我的代码中有任何错误,我想问问您。问题是,我得到了一个解决方案,但它不同于解析解(通过 symbolab 解决)。这是我的代码:
#include <iostream>
#include <vector>
#include <fstream>
using namespace std;
double fx(double x, double y){
return y;
}
double fy(double x, double y){
return x + y;
}
void runge_kutta(double dt, double &x, double &y){
vector<double> k(4), l(4);
#define dt_half dt/2.0
#define dt_six dt/6.0
k.at(0)=fx(x, y);
l.at(0)=fy(x,y);
k.at(1)=fx(x+dt_half*k.at(0), y+dt_half*l.at(0));
l.at(1)=fx(x+dt_half*l.at(0), y+dt_half*l.at(0));
k.at(2)=fx(x+dt_half*k.at(1), y+dt_half*l.at(1));
l.at(2)=fx(x+dt_half*l.at(1), y+dt_half*l.at(1));
k.at(3)=fx(x+dt*k.at(2), y+dt*l.at(2));
l.at(3)=fy(x+dt*k.at(2), y+dt*l.at(2));
x+=dt_six*(k.at(0)+2.0*k.at(1)+2.0*k.at(2)+k.at(3));
y+=dt_six*(l.at(0)+2.0*l.at(1)+2.0*l.at(2)+l.at(3));
}
double t=10.0;
double dt=1e-3;
double x=0.0;
double y=1.0;
ofstream out;
int main(){
out.open("out.txt");
out << "#t, x, y" << endl;
for(double time=0.0;time<t;time+=dt){
cout<<time<<endl;
out << time << ' '
<< x << ' '
<< y << ' '
<< endl;
runge_kutta(dt, x, y);
}
out.close();
}
你能看到任何错误吗?
非常感谢您的帮助!
仔细观察这些线条
l.at(1)=fx(x+dt_half*l.at(0), y+dt_half*l.at(0));
l.at(2)=fx(x+dt_half*l.at(1), y+dt_half*l.at(1));
我认为他们应该使用 fy 而不是 fx。
我已经强调了 RK4 代码中的几个不一致之处
任何处理 y 的东西都应该调用 fy
并且有导数 l
任何处理 x[=22= 的东西] 有导数 k
.
现在我正在努力让我的 Runge-Kutta 发挥作用。实际上我必须解决一个比这个 ODE 系统复杂得多的 ODE 系统,但是 runge-kutta-function 是一样的。我正在解决这个简单的系统,以检查模拟是否收敛到解析解。
首先,如果您发现我的代码中有任何错误,我想问问您。问题是,我得到了一个解决方案,但它不同于解析解(通过 symbolab 解决)。这是我的代码:
#include <iostream>
#include <vector>
#include <fstream>
using namespace std;
double fx(double x, double y){
return y;
}
double fy(double x, double y){
return x + y;
}
void runge_kutta(double dt, double &x, double &y){
vector<double> k(4), l(4);
#define dt_half dt/2.0
#define dt_six dt/6.0
k.at(0)=fx(x, y);
l.at(0)=fy(x,y);
k.at(1)=fx(x+dt_half*k.at(0), y+dt_half*l.at(0));
l.at(1)=fx(x+dt_half*l.at(0), y+dt_half*l.at(0));
k.at(2)=fx(x+dt_half*k.at(1), y+dt_half*l.at(1));
l.at(2)=fx(x+dt_half*l.at(1), y+dt_half*l.at(1));
k.at(3)=fx(x+dt*k.at(2), y+dt*l.at(2));
l.at(3)=fy(x+dt*k.at(2), y+dt*l.at(2));
x+=dt_six*(k.at(0)+2.0*k.at(1)+2.0*k.at(2)+k.at(3));
y+=dt_six*(l.at(0)+2.0*l.at(1)+2.0*l.at(2)+l.at(3));
}
double t=10.0;
double dt=1e-3;
double x=0.0;
double y=1.0;
ofstream out;
int main(){
out.open("out.txt");
out << "#t, x, y" << endl;
for(double time=0.0;time<t;time+=dt){
cout<<time<<endl;
out << time << ' '
<< x << ' '
<< y << ' '
<< endl;
runge_kutta(dt, x, y);
}
out.close();
}
你能看到任何错误吗?
非常感谢您的帮助!
仔细观察这些线条
l.at(1)=fx(x+dt_half*l.at(0), y+dt_half*l.at(0));
l.at(2)=fx(x+dt_half*l.at(1), y+dt_half*l.at(1));
我认为他们应该使用 fy 而不是 fx。
我已经强调了 RK4 代码中的几个不一致之处
任何处理 y 的东西都应该调用 fy
并且有导数 l
任何处理 x[=22= 的东西] 有导数 k
.