我对 C 和 C++ 编程有些陌生,我已经在 R 中实现了它,现在我正在尝试用 C++ 编写它。
我用了std::vector<double>
和 std::vector<std::vector<double>>
并按值传递,因为它一次只需要传递一行来填充 std::vector<std::vector<double>> u(t.size(),vector<double>(n))
其大小基于 h
, k
,以及总运行时间。
我遇到的问题是,当波穿过轴时,输出数据似乎出错了。我无法从逻辑上弄清楚这有什么问题,但我可能漏掉了一些东西,我认为我更有可能滥用了 std::vector
或者有一些我不认识的数据类型冲突。
也许其他人可以看到我看不到的东西,这是我的代码:
#include <iostream>
#include<math.h>
#include<vector>
#include<fstream>
using namespace std;
vector<double> takestep (double h,double k,vector<double> ukm1,vector<double> ukm2){
int n = ukm1.size();
vector<double> uk(n);
uk[0]=0;
uk[n-1]=0;
for(int i = 1; i < (n-1); i++)
{
uk[i] = (pow(k,2)/pow(h,2))*(ukm1[i+1]-2*ukm1[i]+ukm1[i-1])+2*ukm1[i]-ukm2[i];
}
return uk;
}
//-----------------------------------------------------------------
vector<vector<double> > solve1D (double tf, double h, double k, vector<double> ukm1, vector<double> ukm2){
int n = ukm1.size();
vector<double> t((int)tf/k);
for(int i = 0; i<t.size(); i++)
{
t[i] = k*i;
}
vector<vector<double> > u(t.size(),vector<double>(n));
u[0] = ukm1;
u[1] = takestep(k,h,ukm1,ukm2);
for(int i = 2;i<t.size();i++)
{
u[i]= takestep(h,k,u[i-1],u[i-2]);
}
return u;
}
//====================================================================
int main(int argc, char** argv) {
double tf = 12.0;
double k = .005;
double h = .01;
vector<double> x(1.0/h);
for(int i = 1;i<x.size();i++)
{
x[i] = i*h;
}
vector<double> yo(x.size());
yo[0]=0;
yo[x.size()-1]=0;
for(int i = 1;i<yo.size()-1;i++)
{
yo[i] = 0.5*sin(x[i]*M_PI);
}
vector<vector<double> > u = solve1D(tf,k,h,yo,yo);
ofstream myfile;
myfile.open ("Wave1D_output.txt");
for(int i = 0;i<u.size();i++)
{
for(int j = 0;j<yo.size();j++)
{
myfile<< u[i][j]<<"\t";
}
myfile<<"\n";
}
myfile.close();
return 0;
}
我的 R 脚本功能齐全,而且我使用的代码基本相同(除了索引已更改)。主要区别在于 vector 和二维数组(众所周知)在 C++ 中稍微复杂一些。
完成此操作后,我将尝试使用 OpenMP 来完成它,但这(可能)是另一天的问题。
最佳答案
所以,事实证明错误是......我在初始调用中将 h 和 k 向后发送给求解器,我在打印的大小时发现了这一点时间 vector ,它是 1200 而不是 2400。最终保留了原始索引,因为这不是问题所在。如前所述,我添加了一些评论:
#include <iostream>
#include<math.h>
#include<vector>
#include<fstream>
using namespace std;
vector<double> takestep(double h,double k,vector<double> ukm1,vector<double> ukm2){
int n = ukm1.size();//set n to number of elements in of input array
vector<double> uk(n);//new array of same size
uk[0]=0; //
uk[n-1]=0;//set ends to zero
for(int i = 1;i<(n-1);i++){//for i in length of uk
uk[i]=(pow(k,2)/pow(h,2))*(ukm1[i+1]-2*ukm1[i]+ukm1[i-1])+2*ukm1[i]-ukm2[i];//finite difference formula
}
return uk;//returns next time step
}
vector<vector<double> > solve1D (double tf, double h, double k, vector<double> ukm1, vector<double> ukm2){
int n = ukm1.size(); //set n to number of elements in input array
vector<double> t(round(tf/k));//set size of sime sequence (final time/delta time)
for(int i = 0; i<t.size(); i++){
t[i] = k*(i+1); //fill time array with proper sequence
}
vector<vector<double> > u(t.size(),vector<double>(n));//create 2D array to store all the time steps
u[0] = ukm1;//set initial state
u[1] = takestep(k,h,ukm1,ukm2); //set first time step
for(int i = 2;i<t.size();i++){
u[i]= takestep(h,k,u[i-1],u[i-2]);//take all the other time steps
}
return u;//return array with all steps
}
int main(int argc, char** argv) {
double tf = 12.0;//set final time for test
double k = .005;//set time step
double h = .01;//set distance for finite difference
vector<double> x((int)round(1.0/h));//create vector to store x steps
for(int i = 0;i<x.size();i++){
x[i] = i*h; //fill with proper distance steps
}
vector<double> yo(x.size());//create vector to store wave function
for(int i = 1;i<yo.size()-2;i++){
yo[i] = 0.5*sin(x[i]*M_PI); //fill with function values
}
vector<vector<double> > u = solve1D(tf,h,k,yo,yo);//send to solver
关于c++ 1D波函数错误输出(有限差异),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/36094347/