二维热传导方程-有限差分数值算法C++实现

下面的这种二阶偏微分方程即为拉普拉斯方程，在热力学中叫稳态的热传导方程。
$$\frac{{\partial }^{2}U}{\partial x^2}+\frac{{\partial }^{2}U}{\partial y^2}=0$$
二阶差分格式为
$$\frac{U(i+1,j)-2U(i,j)+U(i-1,j)}{d{x}^2}+\frac{U(i,j+1)-2U(i,j)+U(i,j-1)}{d{y}^2}=0$$
当dx和dy为1时，最后化为
$$U(i,j)=\frac{U(i+1,j)+U(i-1,j)+U(i,j-1)+U(i,j+1)}{4}$$
也就是前后左右格子的平均值

下面是一个网格区域，在中间矩形块是温度100的初始值，将偏微分方程写成差分方程，然后执行解算步(sov_step)Ts次以后的格网的温度仿真的可视化结果。

#include <opencv2/opencv.hpp>
#include <opencv2/imgproc/imgproc_c.h>

#include <iostream>
#include<vector>

#include <Eigen/Dense>

void bd_set(Eigen::MatrixXd& m)
{
	int yn = m.rows();
	int xn = m.cols();

	for (int i = 0; i < yn; i++)
	{
		for (int j = 0; j < xn; j++)
		{
			if (i == j)
			{
				m(i, j) = 100.0;
			}
		}
	}
}

void bd_set1(Eigen::MatrixXd& m, int x0, int x1, int y0, int y1)
{
	int yn = m.rows();
	int xn = m.cols();

	for (int i = 0; i < yn; i++)
	{
		for (int j = 0; j < xn; j++)
		{
			if (j >= x0 && j <= x1 &&
				i >= y0 && i <= y1)
			{
				m(i, j) = 100.0;
			}
		}
	}
}

void sov_step(Eigen::MatrixXd& m, int xn, int yn)
{
	for (int i = 1; i < yn; i++)
	{
		for (int j = 1; j <xn; j++)
		{
			if (i == j)
			{
				continue;
			}

			m(i, j) = 0.25 * (m(i - 1, j) + m(i + 1, j) + m(i, j - 1) + m(i, j + 1));
		}
	}
}

void sov_step1(Eigen::MatrixXd& m, int xn, int yn, int x0, int x1, int y0, int y1)
{
	for (int i = 1; i < yn; i++)
	{
		for (int j = 1; j < xn; j++)
		{
			if (j >= x0 && j <= x1 &&
				i >= y0 && i <= y1)
			{
				continue;
			}

			m(i, j) = 0.25 * (m(i - 1, j) + m(i + 1, j) + m(i, j - 1) + m(i, j + 1));
		}
	}
}

int main()
{
	double T = 1.0;

	int nx = 200;
	int ny = 200;


	Eigen::MatrixXd ma = Eigen::MatrixXd::Zero(ny, nx);


	int Ts = 1000;
	int nx_ = nx - 1;
	int ny_ = ny - 1;
	int x0 = nx * 0.3;
	int x1 = nx * 0.7;
	int y0 = ny * 0.3;
	int y1 = ny * 0.7;

	bd_set1(ma, x0,x1,y0,y1);

	for (int k = 0; k < Ts; k++)
	{
		sov_step1(ma, nx_, ny_, x0, x1, y0, y1);
	}

	Mat ma_show = Mat::zeros(nx, ny, CV_8UC3);

	for (int i = 0; i < ny; i++)
	{
		for (int j = 0; j < nx; j++)
		{
			auto& a0 = ma_show.at<cv::Vec3b>(i, j);

			double t0 = ma(i, j) / 100.0;
			double t1 = 1.0 - t0;

			//uchar v = uchar(ma(i,j) * 255.0 / 100.0);
			a0[0] = 80.0 * t0 + 255.0 * t1;
			a0[1] = 240.0 * t0; // 255.0 * a0;
			a0[2] = 255.0 * t0; // 255.0 * a0;
		}
	}

	imshow("image-u", ma_show);

	waitKey();

	destroyAllWindows();

	return 0;
}

仿真结果