tridiag.cpp

追赶法求解三对角矩阵是数值算法中重要的一种,这个利用了求借热传导方程

// How to run the procedure and plot the temperature 
// produce data.dat

// load data.dat
// plot(data) 

/* Tridag for solving Heat conduction problem */

// Solve the tridigonal linear equation

//   1+beta     -1        0          ...   0  |  T^n_0      beta*T_\infty
//   -alpha  1+2 alpha  -alpha       ...   0  |  T^n_1        T^{n-1}_1
//				:			
//				:		   	=
//				:
//     0  ...    0  -alpha  1+2 alpha  -alpha |  T^n_{N-1}    T^{n-1}_{N-1}
//     0  ...    0     0       -1      1+beta |  T^n_N      beta*T_\infty

//   < == >

//     b0       c0        0          ...   0  |    u0           r0
//     a1       b1       c1          ...   0  |    u1           r1
//				:			
//				:		   	=
//				:
//     0  ...    0  a_{N-1}   b_{N-1}  c_{N-1}|  u_{N-1}       r_{N-1}
//     0  ...    0     0        aN         bN |    uN           rN

#include <stdio.h>
#include <iostream.h>
#include <fstream.h>

#define rho 7300 // density
#define Cp 800 // heat capacity
#define K 30 // heat conductivity
#define h 100 // heat exchange coefficient
#define T0 1000 // initial temperature 
#define T_infty 20 // environment temperature 
#define L 0.1 // Fe length


#define dt 1.0 // time step
#define n_step 300 // step number
#define time n_step*dt // Computation will stop at this time from the
		       // initial time 0
#define N 100 // mesh number
#define dx L/N // mesh size

#define alpha (K*dt)/(rho*Cp*dx*dx) 
//#define alpha K*dt/rho/Cp/dx/dx // Wrong,  Who could tell me why?
#define beta (h*dx)/K


void tridag(double a[], double b[], double c[], double r[], double u[], double gam[], int n);

void main(void)
{	int i,step;
	double a[N+1],b[N+1],c[N+1],r[N+1],gam[N+1],T[N+1]; // A[N+1] has no used.
//Initial temperature
	for (i=0;i<=N;i++) T[i]=T0;
//a, b and c will not change with the time iteration
	for (i=1;i<N;i++) {
		a[i]= -alpha;
		b[i]= 1+2*alpha;
		c[i]= -alpha;
	}
	b[0]=1+beta;
	b[N]=1+beta;
	c[0]=-1;
	a[N]=-1;
	
// time evolution
	for (step=0;step<n_step; ++step) {
// form the right hand side(rhs) of equ.
		for (i=1;i<N;i++) r[i]=T[i];
		r[0]=beta*T_infty;
		r[N]=r[0];
// solve this equation		
		tridag(a,b,c,r,T,gam,N+1); // only T changed
	}
// output of data
	for (i=0;i<=N;i++) {
		cout << "T["<< i << "]=" << T[i] <<endl;
	}
	ofstream  os( "data.dat"); 
	for (i=0;i<=N;i++) {
		os << i << " " << T[i] << "\n";
	}

}


void tridag(double a[], double b[], double c[], double r[], double u[], double gam[], int n)
{
        int j;
        double bet; 
                                                                                
        if (b[0] == 0.0) printf("Error 1 in tridag");
        u[0]=r[0]/(bet=b[0]);
        for (j=1;j<n;j++) {
                gam[j]=c[j-1]/bet;
                bet=b[j]-a[j]*gam[j];
                if (bet == 0.0) printf("Error 2 in tridag");
                u[j]=(r[j]-a[j]*u[j-1])/bet;
        }
        for (j=(n-2);j>=0;j--)
                u[j] -= gam[j+1]*u[j+1];
}