z_g_fa.cpp

计算方法程序常微分方程的数值解法课堂讲义

//////////////////////////////////////////////////////////////////////////
//	程序7.5  追赶法求解三对角方程组
//
//////////////////////////////////////////////////////////////////////////

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

#define	MAX_n	100               // 最大阶数

// 解输出
void SulutionOutput( float x[], int n )
{
	int i;

	for( i = 1; i <=n; ++i )
	{
		printf( "\nx[%d] = %f", i , x[i] );
	}
}

// 三对角方程组元素输入
void TriDiagonalMatrixInput( float a[], float b[], float c[], float f[], int n )
{
	int i;

	// 输入第一行
	printf( "Input b[1],c[1],f[1]:" );
	scanf( "%f %f %f", &b[1], &c[1], &f[1] );
	
	// 输入中间各行
	for( i = 2; i < n; ++i )
	{
		printf( "Input a[%d],b[%d],c[%d],f[%d]:", i, i, i, i );
		scanf( "%f %f %f %f", &a[i], &b[i], &c[i], &f[i] );
	}
	
	// 输入最后一行
	printf( "Input a[%d],b[%d],f[%d]:", n, n, n );
	scanf( "%f %f %f", &a[n], &b[n], &f[n] );
}

// 三对角方程组求解 — 追赶法
// 可以在程序中加入当除数为接近0时退出函数的判断
//
void Z_G_method( float a[], float b[], float c[], float f[], int n )
{
	int i;

	//
	c[1] /= b[1];
	for( i = 2; i < n; ++i )
	{
		c[i] /= ( b[i] - a[i] * c[i-1] );
	}

	//
	f[1] /= b[1];
	for( i = 2; i <= n; ++i )
	{
		f[i] = ( f[i] - a[i] * f[i-1] ) / ( b[i] - a[i] * c[i-1] );
	}

	//
	for( i = n-1; i > 0; --i )
	{
		f[i] -= c[i] * f[i+1];
	}
}

void main()
{
	int n;
	float a[MAX_n], b[MAX_n], c[MAX_n], f[MAX_n];
	
	// 输入阶数
	printf( "\nInput n =" );
	scanf( "%d", &n );

	// 输入三对角元素和右端项
	TriDiagonalMatrixInput( a, b, c, f, n );
	
	// 求解
	Z_G_method( a, b, c, f, n );
	
	// 解输出
	SulutionOutput( f, n );
}

/*
	运行实例：

Input n=4
Input b[1],c[1],f[1]:2 1 1
Input a[2],b[2],c[2],f[2]:1 3 1 2
Input a[3],b[3],c[3],f[3]:1 4 1 3
Input a[4],b[4],f[4]:1 5 4

x[1]=0.294118
x[2]=0.411765
x[3]=0.470588
x[4]=0.705882

*/