When trying to plot the solution of the heat PDE I've found some troubles. The solution that I've found is:

Here is the code:

#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;
}

What I'm trying to do is to declare a function that calculates the infinite sum for n odd, and then looping it for every x between 0 and 1. The problem is that I don't know how to declare "result" in order to be the sum of all the terms, because if I declare it outside the for loop it doesn't satisfy the boundary conditions.

(Note that I fixed t=0.025).

CodePudding user response：

According to the equation, you can implement f as:

#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;
}

Since you are using double, so you have to explicitly add a .0 otherwise it may be considered as integer.

The declarations of variable are moved outside the loop in order both to clarify the code and ensure the variable result gets update instead of being overwritten.

EDIT:

You could improve the function f to take the value of t as an input. This also aligns with the equation provided. It would then implements this way:

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;
}

EDIT:

The implementation of the math of the equation could be further simplified:

• a^2 b^2 is same as (ab)^2.
• (-1)^n with n odd is always -1.
• 2*(1-pow((-1),n)) is a replacement for 4.

Plus, from a performance perspective you can avoid recalculation of repeated terms by putting them in a variable and the use it as you need (for instance the n^2 pi^2).