Уменьшение ошибки моей производной функции

Я хотел бы решить формулу

которая описывает комплексную функцию E, продолжающуюся в направлении r (в виде цилиндрической системы координат) и распространяющуюся в направлении z. Из-за невозможности переписать его в производную по времени, я хотел решить это как

используя следующий код:

// MWE_derivation.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"#include "stdafx.h"#define _USE_MATH_DEFINES
#include <cmath>
#include <iostream>
#include <vector>
#include <complex>
#include <fstream>

//Some typedefs
typedef std::complex<double> COMPLEX;
typedef std::vector<COMPLEX> C_VECTOR;
typedef std::vector<C_VECTOR> C_MATRIX;

//Some constants
#define L0 1.2e-6
#define C90 299792458
#define K0 (2 * M_PI / L0)
#define W0 (2 * M_PI*C90 / L0)
#define FWHM0 250e-15
#define W_PULSE 3.25e-6
#define E_VAL 1.60217662e-19
#define ME 9.10938356e-31
#define M0 0.15*ME
//const double rho0 = 0;
#define BETA0 2e-11
#define HBAR 1.0545718e-34
#define N0 3.52
#define N2 4.5e-18
#define ENERGY 2.5e-6
#define POWER 1e7
#define RANGE_TIME 3
#define RANGE_SPACE 3

//Printing function
void print_vector(const C_VECTOR &data, std::string filename)
{
std::ofstream file_out;
file_out.open(filename);
for (unsigned int j = 0; j < data.size(); j++)
file_out << data[j] << ' ';
file_out << '\n';
file_out.close();
}

//Derivative functions

COMPLEX forward_diff_first(const C_VECTOR &x, const int p0, const double dx)
{
return (x[p0 + 1] - x[p0]) / dx;
}

COMPLEX forward_diff_second(const C_VECTOR &x, const int p0, const double dx)
{
return (x[p0 + 2] - x[p0 + 1] - x[p0 + 1] + x[p0]) / (dx*dx);
}

COMPLEX backward_diff_first(const C_VECTOR &x, const int p0, const double dx)
{
return (x[p0] - x[p0 - 1]) / dx;
}

COMPLEX backward_diff_second(const C_VECTOR &x, const int p0, const double dx)
{
return (x[p0 - 2] - x[p0 - 1] - x[p0 - 1] + x[p0]) / (dx*dx);
}

COMPLEX central_diff_first(const C_VECTOR &x, const int p0, const double dx)
{
return (x[p0 + 1] - x[p0 - 1]) / (2 * dx);
}

COMPLEX central_diff_second(const C_VECTOR &x, const int p0, const double dx)
{
return (x[p0 + 2] - x[p0] - x[p0] + x[p0 - 2]) / (4 * dx*dx);
}

C_VECTOR central_diff_f(const C_VECTOR &x, const double dx)
{
C_VECTOR cdf = C_VECTOR(x.size());
cdf[0] = forward_diff_first(x, 0, dx);
cdf[x.size() - 1] = backward_diff_first(x, x.size() - 1, dx);
for (unsigned int i = 1; i < x.size() - 1; i++)
cdf[i] = central_diff_first(x, i, dx);
return cdf;
}

C_VECTOR forward_diff_f(const C_VECTOR &x, const double dx)
{
C_VECTOR cdf = C_VECTOR(x.size());
cdf[x.size() - 1] = backward_diff_first(x, x.size() - 1, dx);
for (unsigned int i = 0; i < x.size() - 1; i++)
cdf[i] = forward_diff_first(x, i, dx);
return cdf;
}

C_VECTOR backward_diff_f(const C_VECTOR &x, const double dx)
{
C_VECTOR cdf = C_VECTOR(x.size());
cdf[0] = forward_diff_first(x, 0, dx);
for (unsigned int i = 1; i < x.size(); i++)
cdf[i] = backward_diff_first(x, i, dx);
return cdf;
}

C_VECTOR central_diff_s(const C_VECTOR &x, const double dx)
{
C_VECTOR cdf = C_VECTOR(x.size());
cdf[x.size() - 1] = backward_diff_second(x, x.size() - 1, dx);
cdf[x.size() - 2] = backward_diff_second(x, x.size() - 2, dx);
cdf[0] = forward_diff_second(x, 0, dx);
cdf[1] = forward_diff_second(x, 1, dx);
for (unsigned int i = 2; i < x.size() - 2; i++)
cdf[i] = central_diff_second(x, i, dx);
return cdf;
}

//Auxiliary functions
bool isnan(const COMPLEX &in)
{
return isnan(in.real()) || isnan(in.imag());
}
template<typename T>
std::vector<T> linspace(T start_in, T end_in, int num_in)
{
double start = static_cast<double>(start_in);
double end = static_cast<double>(end_in);
double num = static_cast<double>(num_in);
double delta = (end - start) / (num - 1);
std::cout << "Temp_data created\n";
std::vector<T> linspaced(num - 1);
for (int i = 0; i < num - 1; ++i)
{
linspaced[i] = start + delta * i;
}
std::cout << "Filled linspaced\n";
linspaced.push_back(end);
std::cout << "linspaced filled finally\n";
return linspaced;
}

C_VECTOR gaussian(const double w, const std::vector<double> r, const double P0)
{
C_VECTOR ret_val(r.size());
for (unsigned int i = 0; i < r.size(); i++)
{
ret_val[i] = P0 * exp(-COMPLEX(2, 0) * r[i] * r[i] / (w*w));
}
return ret_val;
}

double get_abs(COMPLEX number)
{
return number.real()*number.real() + number.imag()*number.imag();
}

//Main function

C_VECTOR deriv_func(const C_VECTOR &in, const double dx)
{
return central_diff_f(in, dx);
}

C_VECTOR calc_steps(const C_VECTOR &in, const std::vector<double> r, const double dz, const int steps, const COMPLEX pref)
{
C_VECTOR calc_val = in;
C_VECTOR add_vec(in.size());
C_VECTOR rho_vals(in.size());
for (unsigned int i = 0; i < in.size(); i++)
rho_vals[i] = 1e20;
for (int i = 0; i < steps; i++)
{
print_vector(calc_val, "new_vec_" + std::to_string(i) + ".txt");
C_VECTOR deriv = deriv_func(calc_val, abs(r[0] - r[1]));
C_VECTOR deriv2 = central_diff_s(calc_val, abs(r[0] - r[1]));
for (unsigned int j = 1; j < in.size(); j++)
{
if (isnan(deriv[j]))
{
std::cout << "isnan value in round " << i << " detected.\n ";
return calc_val;
}
}

for (unsigned int j = 0; j < in.size(); j++)
{
double r_val = r[j] + 1e-9;
add_vec[j] = pref * dz*(deriv[j] / r_val + deriv2[j])* calc_val[j] + COMPLEX(0, 1)*(K0)*N2 / N0*get_abs(calc_val[j])*calc_val[j] - BETA0 / 2 * calc_val[j] * get_abs(calc_val[j]) - 2*M_PI*COMPLEX(0, 1)*K0*E_VAL*E_VAL/(N0*N0*M0*W0*W0)*rho_vals[j]*calc_val[j];
calc_val[j] += add_vec[j];
}
//calc_val[0] = pref * dz * (deriv[1] / r[1] + deriv2[0])*calc_val[0];}
return calc_val;
}

int main()
{
double w = 3.25e-6;
int num_vals = 1000;
std::vector<double> r_val = linspace<double>(0.0, 3 * w, num_vals);
C_VECTOR gauss = gaussian(w, r_val, 1e-2);
print_vector(gauss, "gauss.txt");
std::cout << "pref is " << COMPLEX(0, 1 / (2 * K0)) << '\n';
C_VECTOR new_vec = calc_steps(gauss, r_val, 1e-7, 1000, COMPLEX(0, 1)/(2*K0));
print_vector(new_vec, "new_vec.txt");
return 0;
}

К сожалению, у меня проблема в том, что мои производные функции добавляют артефакты к функции, что приводит к перерегулированию и чрезвычайно большим числам на границах. Это можно увидеть в приведенном выше коде с текущими значениями и уже в итерации 2 (в файле new_vec_2.txt, на левой границе). Полученная функция больше не является гладкой. Добавление большего количества баллов только задерживает эту проблему, но не решает ее. Что я могу сделать, чтобы улучшить код и предотвратить перерегулирование? Или алгоритм неисправен, и мне нужно разработать другой?