C++实现前向欧拉法Forward Euler解决偏微分方程

1. 数学原理

以简单的热方程为例，其中D为常数：

根据导数的定义式有：

如果我们将x理解成delta_x*n，则u(x,t)可以看成只是t的函数，n作为参数来控制x，即定义：

这样我们就可以把这个偏微分方程看做是一个常微分方程

对于常微分方程的C++实现，可以先看看https://blog.csdn.net/weixin_39374967/article/details/105333880

二、C++实现

如下代码中，新建了一个叫solutionDirichletProblem.csv的文件记录偏微分方程的计算结果；

初始条件是

#define _USE_MATH_DEFINES
#include <cmath>
#include <iostream>
#include <fstream>
#include <vector>

// 设置左右边界条件
double DirichletBCLeft(double t)
{
	return 0.0;
}
double DirichletBCRight(double t)
{
	return 0.0;
}

// 新建一个文件，记录时间和函数值
void AddCurrentStateToOpenStream(double currentTime, const std::vector<double>& currentFunctionValue)
{
	std::ofstream fileOutput("solutionDirichletProblem.csv", std::ofstream::app);
	fileOutput << currentTime << " , ";

	for (int i = 0; i < currentFunctionValue.size() - 1; ++i)
	{
		fileOutput << currentFunctionValue[i] << " , ";
	}

	fileOutput << currentFunctionValue.back() << std::endl;

	fileOutput.close();
}

int main()
{
	const double finalTime = 0.1;
	const int numberOfTimeSteps = 10000;
	const double deltaT = finalTime / numberOfTimeSteps;

	// x的取值范围是[0,1]
	const double domainSize = 1.0;
	const int numberOfSpatialElements = 51;
	const double deltaX = domainSize / (numberOfSpatialElements - 1);

	// 初始条件里，开始结尾都是0，中间部分是sin(pi*x)
	std::vector<double> initialFunctionValue(numberOfSpatialElements, 0.0);
	initialFunctionValue[0] = DirichletBCLeft(0.0);
	initialFunctionValue.back() = DirichletBCRight(0.0);
	for (int i = 1; i < numberOfSpatialElements - 1; ++i)
	{
		initialFunctionValue[i] = sin(M_PI*i*deltaX);
	}

	// initial是t=0时的x序列，previous是t时的x序列，current是t+deltaT的序列
	std::vector<double> previousFunctionValue = initialFunctionValue;
	std::vector<double> currentFunctionValue = previousFunctionValue;

	double currentTime = 0.0;

	double outputTimesInterval = finalTime / 10;
	double nextOutputTime = outputTimesInterval;
	AddCurrentStateToOpenStream(currentTime, previousFunctionValue);

	while (currentTime < finalTime)
	{
		previousFunctionValue = currentFunctionValue;

		for (int i = 1; i < numberOfSpatialElements - 1; ++i)
		{
			currentFunctionValue[i] += deltaT / deltaX / deltaX * (previousFunctionValue[i + 1] + previousFunctionValue[i - 1] - 2 * previousFunctionValue[i]);
		}

		// 更新开头结尾边界条件
		currentFunctionValue[0] = DirichletBCLeft(currentTime);
		currentFunctionValue.back() = DirichletBCRight(currentTime);

		if (currentTime >= nextOutputTime)
		{
			AddCurrentStateToOpenStream(currentTime, currentFunctionValue);
			nextOutputTime += outputTimesInterval;
		}

		currentTime += deltaT;
	}

	return 0;
}