热方程解

Question

当我试图绘制出热PDE的解时，我发现了一些麻烦。我找到了是的解决方案

以下是代码：

#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define N 10000
 
double f(double x);
double X[N];
double Y[N];
 
int main(){
    int i;
    double b=43351/94400;
    double dx=0.0001;
    X[0]=0;
    Y[0]=b;
    for (i=1; i<N; i++){
    X[i]=X[i-1]+dx;
    Y[i]=f(X[i]);
    }
    FILE* output;
    
    output = fopen("dades.txt", "w");
   
    
    fprintf(output, "x Posició Temperatura\n");
    
    for (i = 0; i < N; i++){
        fprintf(output, "%lf %lf %lf\n", i*dx, X[i], Y[i]);
    }
    fclose(output);
    return 0;
}
 
double f(double x){
    
    int n;
    double b=43351/94400;
    for (n=1; n<N; n+=2){
     
        double pi=3.14159265358979323846;
        double t0=0.025;
        double result=b;
        result+=2*(1-pow((-1),n))/(pi*n)*(1-exp(-pow(n,2)*pow(pi,2)*pow(t0,2)))/(pow(n,2)*pow(pi,2))*sin(n*pi*x);
        }
    return result;
}

我要做的是声明一个函数，它计算n个奇数的无穷和，然后把它循环到0到1之间的每一个x。问题是，我不知道如何声明“结果”来作为所有项的和，因为如果我在for循环之外声明它，它不满足边界条件。

(请注意，我修复了t=0.025)。

Answer 1

根据该公式，您可以将f实现为：

#define M_PI 3.14159265358979323846;

double f(double x)
{
    int n;
    double result=43351.0/94400.0;
    double t0=0.025;

    for (n=1; n<N; n+=2){
        result+=2*(1-pow((-1),n))/(M_PI*n)*(1-exp(-pow(n,2)*pow(M_PI,2)*pow(t0,2)))/(pow(n,2)*pow(M_PI,2))*sin(n*M_PI*x);
    }
    return result;
}

因为您使用的是double，所以您必须显式地添加一个.0，否则它可能被认为是整数。

变量的声明被移出循环之外，以便澄清代码并确保变量result得到更新而不是被覆盖。

编辑：

您可以改进函数f，将t的值作为输入。这也符合所提供的方程式。然后，它将以这种方式实现：

double f(double x, double t)
{
    int n;
    double result=43351.0/94400.0;

    for (n=1; n<N; n+=2){
        result+=2*(1-pow((-1),n))/(M_PI*n)*(1-exp(-pow(n,2)*pow(M_PI,2)*pow(t,2)))/(pow(n,2)*pow(M_PI,2))*sin(n*M_PI*x);
    }
    return result;
}

编辑：

方程式的数学实现可以进一步简化：

• a^2 b^2和(ab)^2是一样的。
• N个奇数的(-1)^n总是-1。
• 2*(1-pow((-1),n))是4的替代品。

此外，从性能的角度来看，您可以通过将重复的术语放入变量并根据需要使用它们(例如，n^2 pi^2)来避免重复的术语的重新计算。