Недавно я снова взял программирование для школы, и у меня возникла проблема с моим кодом, которая включает в себя пакет быстрого преобразования Фурье FFTW.
В моем коде я начинаю с начальной функции (в этом случае u (x, t = 0) = sin (x) + sin (3 * x) и буду использовать RK4, чтобы попытаться решить U_t уравнения теплопроводности Для тех, кто имеет опыт работы с пакетом FFTW, я отправляю массив размера N и преобразовываю его, беря две производные, а затем выполняю обратное преобразование, чтобы найти U_xx уравнения теплопроводности.
#include <iostream>
#include <cmath>
#include <fstream>
#include <fftw3.h>
using namespace std;
const double Pi= 3.141592653589;
const int N = 2048;
void Derivative(double& num1, double& num2, int J){
//takes two derivatives
num1 = (-1)*J*J*num1;
num2 = (-1)*J*J*num2;
}
double Function(double X){ // Initial function
double value = sin(X) + sin(3*X);
return value;
}
void FFT(double *in, double *Result){
double *input, *outI;
input = new double [N];
outI = new double [N];
for(int i=0; i<N; i++){
input[i] = in[i];
}
//Declaring transformed matricies;
fftw_complex *out;
out= (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*((N/2)+1));
//Plans for execution
fftw_plan p;
p= fftw_plan_dft_r2c_1d(N,input, out, FFTW_ESTIMATE);
fftw_execute(p);
//Derivative
for(int i=0; i<(N/2+1);i++){
Derivative(out[i][0],out[i][1],i);
}
fftw_plan pI;
//Execution of Inverse FFT
pI = fftw_plan_dft_c2r_1d(N,out,outI,FFTW_PRESERVE_INPUT);
fftw_execute(pI);
for(int i=0;i<N; i++){//Dividing by the size of the array
Result[i] = (1.0/N)*outI[i];
}
fftw_free(outI); fftw_free(out);
fftw_free(input);
fftw_destroy_plan(pI); fftw_destroy_plan(p);
}
int main()
{
//Creating and delcaring function variables
double *initial,*w1,*w2,*w3,*w4,*y,*holder1,*holder2, *holder3;
double dx= 2*Pi/N, x=0, t=0;
double h = pow(dx,2), tmax= 100*h;
//creating arrays for the RK4 procedure
initial = new double [N];
w1 = new double [N]; w2 = new double [N]; y= new double [N];
w3 = new double [N]; w4 = new double [N]; holder1 = new double [N];
holder2 = new double [N]; holder3 = new double [N];
double *resultArray=new double[N];
int j =0;
//Output files
ofstream sendtofileINITIAL("HeatdataINITIAL.dat");
ofstream sendtofile5("Heatdata5.dat");
ofstream sendtofile25("Heatdata25.dat");
ofstream sendtofile85("Heatdata85.dat");
ofstream sendtofileFINAL("HeatdataFINAL.dat");
//Initial data
for(int i=0; i<N; i++){
y[i] = Function(x);
sendtofileINITIAL << i*2*Pi/N <<" "<< y[i] << endl;
x+=dx;
}
//RK4
// y[i+1] = y[i] + (h/2)*(w1 + 2*w2 + 2*w3 + w4)
// where the w1,w2,w3,w4 is the standard RK4 calculations
// w1= f(y[i]), w2= f(y[i]+(h/2)*w1[i])
// w3= f(y[i] + (h/2)*w2[i]) and w4 = f(y[i]+h*w3[i])
while(t<=tmax){
FFT(y,resultArray);
for(int i=0; i<N;i++){ //Calculating w1
w1[i] = resultArray[i];
}
for(int i=0; i<N; i++){
holder1[i] = y[i] + (h/2)*w1[i];}
FFT(holder1,resultArray);
for(int i=0; i<N; i++){ //Calculating w2
w2[i] = resultArray[i];
}
for(int i=0; i<N; i++){
holder2[i] = y[i] + (h/2)*w2[i];
}
FFT(holder2,resultArray);
for(int i=0; i<N; i++){ //Calculating w3
w3[i] = resultArray[i];
}
for(int i=0; i<N; i++){
holder3[i]= y[i] + h*w3[i];
}
FFT(holder3,resultArray);
for(int i=0; i<N; i++){ //Calculating w4
w4[i] = resultArray[i];
}
for(int i=0; i<N; i++){
y[i] += (h/6.0)*(w1[i] + 2*w2[i] + 2*w3[i] +w4[i]);
//Outputing data at certain time periods
if(j==5)
sendtofile5 << i*2*Pi/N <<" "<< y[i] << endl;
if(j==25)
sendtofile25 << i*2*Pi/N <<" "<< y[i] << endl;
if(j==85)
sendtofile85 << i*2*Pi/N <<" "<< y[i] << endl;
if(j==100)
sendtofileFINAL << i*2*Pi/N <<" "<< y[i] << endl;
}
j++;// counter
t+=h;// time counter
}
sendtofileINITIAL.close();
sendtofile5.close();
sendtofile25.close();
sendtofile85.close();
sendtofileFINAL.close();
return 0;
}
Когда я строю график y [i] после нескольких итераций t, значения взрываются довольно быстро и чередуются между положительным и отрицательным.
Если у кого-то есть какие-либо предложения по синтаксическим проблемам, которые могут вызывать эти странные результаты, или если я обнаружу ошибки в моих математических методах, я был бы очень признателен за любую информацию. Спасибо.