一道《微分方程数值解法》作业题

/*****************************************************************************
FileName: chafen.c
About: 解决微分方程数值解的一道题目 
Author: Alf
Email: naihe2010@126.com
Url: http://alf.ik8.com
*****************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <error.h>
#include <sys/types.h>
#include <unistd.h>
#include <fcntl.h>

/* 常量定义，可以更改进行试验
   其实t可以设得尽量小一些，但为了尽快出结果，我就先设成这么大 */
#define b 1
#define r 0.5
#define t 0.1
#define h t/r


/* 定义t向及n向节点数 */
int T = (int)(1 / t + 1), H = (int)(1 / h + 1);


/* 定义u(x,0)={}等式 */
#define u(x) ((x<=0.5)? (x): (1-x))

/* 定义显示输出并写入文件过程 */
int wriout(FILE *fp, char *output, double *u){  
    int i;

    if ( (fp = fopen(output, "a+")) == NULL ){
	printf("open file: %s, error!\n", output);
	exit(1);
    }
    /* 我没用过Matlab，不知道输入文件应该是什么格式，只能暂时写成一个数字、一个制表符的循环 */
    for(i=0; i<H; ++i){
	fprintf ( fp, "%f\t", u[i] );
    }
    
    fclose(fp);
    return 0;
}

int main(int argc, char *argv[])
{
    int n, j;
    FILE *fp1, *fp2;
    double u1[H], u2[H];
    
    /* 当n=0时，使用u(x,0)等式 */
    for (j=0; j<H; ++j) {
       	u1[j] = u(j*t);
	u2[j] = u(j*t);
    }
    /* 分别输出到文件中 */
    wriout(fp1, "out1", u1);
    wriout(fp1, "out2", u2);

    /* 当n>=1时，使用差分格式 */
    for(n=1; n<T; ++n){

	/* 方程中第三个等式。
	   其实对于较小的t及h时，这两步完全可以不要。
	   因为后面的叠代过程显然都没有对它们造成影响 */
	u1[0] = u1[H-1] = 0;
	u2[0] = u2[H-1] = 0;

	for(j=1; j<H-1; ++j){
	    /* 第一个差分格式 */
	    u1[j] = u1[j] \
		+ t * (u1[j+1] - 2 * u1[j] + u1[j-1]) / ( h * h ) \
		- t * b * (u1[j+1] - u1[j-1]) / ( 2 * h );
	    /* 另一个差分格式 */
	    u2[j] = u2[j] \
		+ t * (u2[j+1] - 2 * u2[j] + u2[j-1]) / ( h * h ) \
		- t * b * (u2[j] - u2[j-1]) / h;
	}

	/* 因为下一次叠代将失去这组值，所以，只能先保存到文件中，
	   同理，每次n++，都得把前一组先保存 */
	wriout(fp1, "out1", u1);
	wriout(fp2, "out2", u2);
    }
    return 0;
}