c++ 一维波函数坏输出(有限差分)

Question

我对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 脚本功能齐全，我使用的代码几乎相同（除了更改索引）。主要区别在于向量和 2D 数组（如您所知）在C++中稍微复杂一些。

完成此操作后，我将尝试使用 OpenMP 执行此操作，但这（可能）是另一天的另一个问题。

Answer 1

所以，事实证明这个错误是......我在初始调用中向后发送 h 和 k 到求解器，当我打印时间向量的大小时我发现了这一点，它是 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