在C语言中对矩阵应用梯形积分法的实现疑问

问题背景 #

我有一个名为one_elect的2x2矩阵，它是变量r的函数。作为C语言入门学习者，不清楚如何对该矩阵应用梯形积分法，也不知道如何针对变量r循环处理矩阵。以下是我编写的代码和梯形积分法示例代码：

原代码（存在问题） #

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define Z 1
#define BASISFUN 2
#define UPPER 100 // 积分上限
#define LOWER 0   // 积分下限
#define SUBINT 50 // 积分区间数

void makingmatrix(double a[][BASISFUN], double d[][BASISFUN], double alpha_coeff[BASISFUN], double d_coeff[BASISFUN]){
    
    for (int i = 0; i < BASISFUN; i++){
        for (int j = 0; j < BASISFUN; j++){
            a[i][j] = alpha_coeff[j];
            d[i][j] = d_coeff[j];
        }
    }

    printf("\nAlpha 2x2 Coefficient Matrix is:\n");
    for (int i = 0; i < BASISFUN; i++){
        for (int j = 0; j < BASISFUN; j++){
            printf("%.4f\t", a[i][j]);
        }
        printf("\n");
    }

   printf("\nD 2x2 Coefficient Matrix is:\n");
    for (int i = 0; i < BASISFUN; i++){
        for (int j = 0; j < BASISFUN; j++){
            printf("%.4f\t", d[i][j]);
        }
        printf("\n");    
    }
    printf("\n------------------------------------------------------------------\n");
};

void one_ele(double a[][BASISFUN], double d[][BASISFUN], double r, double one_electron[][BASISFUN]){

    double i_step = (UPPER-LOWER)/SUBINT;

    for (int i = 0; i < BASISFUN; i++){
            for (int j = 0; j < BASISFUN; j++){
                
                if (i != j){
                    one_electron[i][j] = (0.5) * (-2 * a[i][j] * r * exp(-a[i][j] * pow(r, 2))) * (-2 * a[j][i] * r * exp(-a[j][i] * pow(r, 2))) + (Z / r);
                } 
                
                else{
                    one_electron[i][i] = (0.5) * (-2 * a[i][i] * r * exp(-a[i][i] * pow(r, 2))) * (-2 * a[i][i] * r * exp(-a[i][i] * pow(r, 2))) + (Z / r);
                }
            }
    }
 printf("\nOne Electron Hamiltonian Elements are\n");
    for (int i = 0; i < BASISFUN; i++){
        for (int j = 0; j < BASISFUN; j++){
            printf("%.4f\t", one_electron[i][j]);
        }
        printf("\n");
    }

    printf("\n------------------------------------------------------------------\n");
}

int main(){
    double A1 = 0.532149, A2 = 4.097728, D1 = 0.82559, D2 = 0.28317;
    double r;
    
    double alpha_coeff[BASISFUN] = {A1, A2};
    double d_coeff[BASISFUN] = {D1, D2};
    double a[BASISFUN][BASISFUN], d[BASISFUN][BASISFUN];
    
    makingmatrix(a, d, alpha_coeff, d_coeff);
    one_ele(a, d, r, one_electron);
    
    return 0;
}

梯形积分示例代码（存在问题） #

#include <stdio.h>
#include <math.h>

#define BASISFUN 2 
#define Z 1.0
#define UPPER 100 // 积分上限
#define LOWER 0   // 积分下限
#define SUBINT 50 // 积分区间数

double f(double r) {
    // 定义函数
    return (0.5) * (-2 * a[i][j] * r * exp(-a[i][j] * pow(r, 2))) * (-2 * a[i][j] * r * exp(-a[i][j] * pow(r, 2))) + (Z / r);
}

int main() {
    double stepSize = (UPPER - LOWER) / SUBINT;
    double integral = f(a) + f(b);

    for (int i = 1; i < SUBINT; i++) {
        double r = a + i * stepSize;
        integral += 2 * f(r);
    }

    integral *= stepSize / 2.0;

    printf("Required value of integration is: %.3f\n", integral);
    return 0;
}

修正后的完整实现 #

核心思路：对矩阵的每个元素单独执行梯形积分，因为每个元素都是r的函数。同时修复原代码中的变量未声明、奇点（r=0时除以0）等问题。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define Z 1
#define BASISFUN 2
#define UPPER 100.0   // 积分上限
#define LOWER 1e-6    // 修改为极小值避免r=0的奇点
#define SUBINT 500    // 增加区间数提升精度

// 生成系数矩阵
void makingmatrix(double a[][BASISFUN], double d[][BASISFUN], double alpha_coeff[BASISFUN], double d_coeff[BASISFUN]){
    for (int i = 0; i < BASISFUN; i++){
        for (int j = 0; j < BASISFUN; j++){
            a[i][j] = alpha_coeff[j];
            d[i][j] = d_coeff[j];
        }
    }

    printf("\nAlpha 2x2 Coefficient Matrix:\n");
    for (int i = 0; i < BASISFUN; i++){
        for (int j = 0; j < BASISFUN; j++){
            printf("%.6f\t", a[i][j]);
        }
        printf("\n");
    }

    printf("\nD 2x2 Coefficient Matrix:\n");
    for (int i = 0; i < BASISFUN; i++){
        for (int j = 0; j < BASISFUN; j++){
            printf("%.6f\t", d[i][j]);
        }
        printf("\n");    
    }
    printf("\n---------------------------------------------------------\n");
}

// 计算指定(i,j)位置的矩阵元素在r处的值
double compute_matrix_element(double a[][BASISFUN], int i, int j, double r){
    double term1, term2;
    if (i != j) {
        term1 = (-2 * a[i][j] * r * exp(-a[i][j] * pow(r, 2)));
        term2 = (-2 * a[j][i] * r * exp(-a[j][i] * pow(r, 2)));
    } else {
        term1 = (-2 * a[i][i] * r * exp(-a[i][i] * pow(r, 2)));
        term2 = term1; // 对角元素两个项相同
    }
    return 0.5 * term1 * term2 + (Z / r);
}

// 对整个矩阵应用梯形积分，结果存入integrated_matrix
void trapz_integrate_matrix(double a[][BASISFUN], double integrated_matrix[][BASISFUN]){
    double step_size = (UPPER - LOWER) / SUBINT;

    // 初始化积分矩阵
    for (int i = 0; i < BASISFUN; i++){
        for (int j = 0; j < BASISFUN; j++){
            // 梯形积分第一步：上下限的函数值之和
            integrated_matrix[i][j] = compute_matrix_element(a, i, j, LOWER) + compute_matrix_element(a, i, j, UPPER);
        }
    }

    // 遍历中间区间
    for (int k = 1; k < SUBINT; k++){
        double r = LOWER + k * step_size;
        for (int i = 0; i < BASISFUN; i++){
            for (int j = 0; j < BASISFUN; j++){
                integrated_matrix[i][j] += 2 * compute_matrix_element(a, i, j, r);
            }
        }
    }

    // 乘以步长的一半完成积分
    for (int i = 0; i < BASISFUN; i++){
        for (int j = 0; j < BASISFUN; j++){
            integrated_matrix[i][j] *= step_size / 2.0;
        }
    }
}

// 打印矩阵
void print_matrix(double mat[][BASISFUN], const char* name){
    printf("\n%s:\n", name);
    for (int i = 0; i < BASISFUN; i++){
        for (int j = 0; j < BASISFUN; j++){
            printf("%.6f\t", mat[i][j]);
        }
        printf("\n");
    }
    printf("\n");
}

int main(){
    double A1 = 0.532149, A2 = 4.097728, D1 = 0.82559, D2 = 0.28317;
    
    double alpha_coeff[BASISFUN] = {A1, A2};
    double d_coeff[BASISFUN] = {D1, D2};
    double a[BASISFUN][BASISFUN], d[BASISFUN][BASISFUN];
    double integrated_one_elect[BASISFUN][BASISFUN]; // 积分后的矩阵

    makingmatrix(a, d, alpha_coeff, d_coeff);
    trapz_integrate_matrix(a, integrated_one_elect);
    print_matrix(integrated_one_elect, "Integrated One-Electron Hamiltonian Matrix");

    return 0;
}

关键修改说明 #

• 修复奇点问题：将积分下限从0改为1e-6，避免Z/r在r=0时出现除以0的错误。
• 拆分功能函数：
  • compute_matrix_element：单独计算矩阵某个位置(i,j)在特定r下的值，便于在积分循环中重复调用。
  • trapz_integrate_matrix：对矩阵每个元素应用梯形积分，外层循环遍历r的子区间，内层循环遍历矩阵的每个元素。
• 修正变量错误：原代码中main函数未声明one_electron变量，示例代码中a、b、i、j未定义，均已修复。
• 提升精度：将子区间数从50改为500，平衡计算速度和积分精度。

内容的提问来源于stack exchange，提问作者Moe El