长方体3维热传导偏微分方程C++算法实现，其中要求长方体表面和外界采用辐射换热

长方体的热传导偏微分方程可以写成：

ρc ∂T/∂t = k (∂²T/∂x² + ∂²T/∂y² + ∂²T/∂z²)

其中，ρ是密度，c是比热容，k是导热系数，T是温度，t是时间，x、y、z是空间坐标。

为了实现表面和外界的辐射换热，可以采用下面的算法：

1. 初始化长方体内部每个点的温度，以及长方体表面和外界的温度。

2. 按照时间步长，进行迭代计算。每个时间步长内，按照以下步骤进行：

2.1 计算长方体内部每个点的温度变化量，根据上述偏微分方程计算。其中，需要注意边界条件，即长方体表面和外界的温度对内部温度的影响。

2.2 根据计算出的温度变化量，更新长方体内部每个点的温度。

2.3 更新长方体表面和外界的温度，根据辐射换热的公式计算。具体公式可以根据实际情况进行选择。常用的有斯特凡-玻尔兹曼定律、维恩位移定律等。

3. 根据需要输出长方体内部每个点的温度分布图，以及长方体表面和外界的温度变化图。

下面是一个简单的C++实现代码，供参考：

#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;

const int N = 10; //长方体边长
const double rho = 1.0; //密度
const double c = 1.0; //比热容
const double k = 1.0; //导热系数
const double dt = 0.01; //时间步长
const double dx = 0.1; //空间步长

double T[N+2][N+2][N+2]; //温度
double T_surf[N+2][N+2]; //表面温度

//初始化
void init() {
    for (int i = 0; i < N+2; i++) {
        for (int j = 0; j < N+2; j++) {
            for (int k = 0; k < N+2; k++) {
                T[i][j][k] = 0.0;
            }
        }
    }
    for (int i = 0; i < N+2; i++) {
        for (int j = 0; j < N+2; j++) {
            T_surf[i][j] = 20.0;
        }
    }
}

//计算温度变化量
double calc_delta_T(int i, int j, int k) {
    double dT = 0.0;
    dT += k * (T[i+1][j][k] - 2.0*T[i][j][k] + T[i-1][j][k]) / pow(dx, 2);
    dT += k * (T[i][j+1][k] - 2.0*T[i][j][k] + T[i][j-1][k]) / pow(dx, 2);
    dT += k * (T[i][j][k+1] - 2.0*T[i][j][k] + T[i][j][k-1]) / pow(dx, 2);
    dT *= dt / (rho * c);
    return dT;
}

//更新温度
void update_T() {
    for (int i = 1; i < N+1; i++) {
        for (int j = 1; j < N+1; j++) {
            for (int k = 1; k < N+1; k++) {
                T[i][j][k] += calc_delta_T(i, j, k);
            }
        }
    }
}

//计算表面温度变化量
double calc_delta_T_surf(int i, int j) {
    double dT_surf = 0.0;
    dT_surf += 4.0 * k * (T_surf[i+1][j] - T_surf[i][j]) / pow(dx, 2);
    dT_surf += 4.0 * k * (T_surf[i-1][j] - T_surf[i][j]) / pow(dx, 2);
    dT_surf += 4.0 * k * (T_surf[i][j+1] - T_surf[i][j]) / pow(dx, 2);
    dT_surf += 4.0 * k * (T_surf[i][j-1] - T_surf[i][j]) / pow(dx, 2);
    dT_surf *= dt / (rho * c);
    return dT_surf;
}

//更新表面温度
void update_T_surf() {
    for (int i = 1; i < N+1; i++) {
        for (int j = 1; j < N+1; j++) {
            T_surf[i][j] += calc_delta_T_surf(i, j);
        }
    }
}

//输出温度分布图
void output_T(int t) {
    char filename[20];
    sprintf(filename, "T%d.dat", t);
    ofstream fout(filename);
    for (int i = 1; i < N+1; i++) {
        for (int j = 1; j < N+1; j++) {
            for (int k = 1; k < N+1; k++) {
                fout << T[i][j][k] << " ";
            }
            fout << endl;
        }
        fout << endl;
    }
    fout.close();
}

//输出表面温度变化图
void output_T_surf(int t) {
    char filename[20];
    sprintf(filename, "T_surf%d.dat", t);
    ofstream fout(filename);
    for (int i = 1; i < N+1; i++) {
        for (int j = 1; j < N+1; j++) {
            fout << T_surf[i][j] << " ";
        }
        fout << endl;
    }
    fout.close();
}

int main() {
    init();
    for (int t = 0; t < 100; t++) {
        update_T();
        update_T_surf();
        output_T(t);
        output_T_surf(t);
    }
    return 0;
}