Next: 4 C & GLSC Up: 3 C & fplot Previous: 3.2 陽解法 heat1d-e.c
/*
* heat1d-i.c -- 陰的スキームで熱方程式を解く
*
* 菊地・山本, 微分方程式と計算機演習, 山海堂 (1991) 第11章
* p237 のプログラムを修正・拡張したもの
*
C ** IMPLICIT FINITE DIFFERENCE METHOD **
C ** FOR HEAT EQUATION **
C nfunc : 関数の番号（5種類の関数を用意している）
C theta : スキームのパラメーター
C N : 分割の数
C t : 時刻
C tau : 時間の刻み幅
C Tmax : 時刻の上限
C h : 空間の刻み幅
C AL,AD,AU : 係数行列
C u : DIMENSION FOR UNKNOWN FUNCTION
*
* このプログラムは
* http://nalab.mind.meiji.ac.jp/%7Emk/program/
* から入手可能
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fplot.h>
/* 円周率 --- main() で pi = 4.0 * atan(1.0); と代入 */
double pi;
void trilu(int, double *, double *, double *);
void trisol(int, double *, double *, double *, double *);
double f(double, int);
void print100(int, double, double *);
void axis(double, double, double, double);
void fsymbol(double, double, char *);
void print_parameters(double, double, double, double,
int, double, int,
double, double, double, double);
void save_picture(void);
int main()
{
int N, nfunc, i, n, nMax;
double a, b, theta, h, Tmax, tau, c1, c2, c3, c4, lambda, t;
double *AL, *AD, *AU, *ff, *u;
double xleft, ybottom, xright, ytop;
/* グラフを書き換える時間間隔 */
double dt;
/* 何ステップおきにグラフを書き換えるか */
int skip;
/* 以前描いたグラフを消すか？ (0=No, 1=Yes) */
int erase_always = 0;
/* 数値を表示するか? (0=No, 1=Yes) */
int print_numerical_data = 0;
/* 円周率 */
pi = 4.0 * atan(1.0);
/* 区間の左端、右端 */
a = 0.0;
b = 1.0;
/* 入力 */
printf("入力して下さい : nfunc(1..5)=");
scanf("%d", &nfunc);
printf("入力して下さい : θ=");
scanf("%lf", &theta);
printf("入力して下さい : N=");
scanf("%d", &N);
printf("入力して下さい : λ=");
scanf("%lf", &lambda);
/* ベクトルの準備 */
if ((AL = malloc(N * sizeof(double))) == NULL) {
fprintf(stderr, "cannot allocate AL[]\n");
exit(1);
}
if ((AD = malloc(N * sizeof(double))) == NULL) {
fprintf(stderr, "cannot allocate AD[]\n");
exit(1);
}
if ((AU = malloc(N * sizeof(double))) == NULL) {
fprintf(stderr, "cannot allocate AU[]\n");
exit(1);
}
if ((ff = malloc(N * sizeof(double))) == NULL) {
fprintf(stderr, "cannot allocate ff[]\n");
exit(1);
}
if ((u = malloc((N + 1) * sizeof(double))) == NULL) {
fprintf(stderr, "cannot allocate u[]\n");
exit(1);
}
/* */
h = (b - a) / N;
tau = lambda * h * h;
printf(" 時間の刻み幅τ= %g になりました。\n", tau);
/* 入力 */
printf("入力して下さい : 最終時刻 Tmax=");
scanf("%lf", &Tmax);
printf("入力して下さい : グラフ書き換え時間間隔 Δt=");
scanf("%lf", &dt);
skip = (int) rint(dt / tau);
/* skip が 0 になっていないかチェック */
if (skip == 0) {
fprintf(stderr, "Δtが小さすぎるので、Δt=τ=%g にします。\n", tau);
skip = 1;
dt = tau;
}
/* 係数行列の準備 */
c1 = - theta * lambda;
c2 = 1.0 + 2.0 * theta * lambda;
c3 = 1.0 - 2.0 * (1.0 - theta) * lambda;
c4 = (1.0 - theta) * lambda;
for (i = 1; i < N; i++) {
AL[i] = c1;
AD[i] = c2;
AU[i] = c1;
}
/* 初期値 */
for (i = 0; i <= N; i++)
u[i] = f(a + i * h, nfunc);
/* 出力（t=0） */
/* 出力（ここでは画面に数値を表示すること）は FORTRAN と C で
かなり異なるので、直接の翻訳は出来ない。ごたごたするので、
print100 という関数にまとめた */
t = 0.0;
print100(N, t, u);
/* ***************** グラフィックスの準備 ***************** */
/* グラフィックス画面の準備 -- まず画面の範囲を決めてから */
xleft = a - (b - a) / 10;
xright = b + (b - a) / 10;
ybottom = - 0.1;
ytop = 1.1;
openpl();
fspace2(xleft, ybottom, xright, ytop);
/* パラメーター、入力データ等の表示 */
print_parameters(xleft, ybottom, xright, ytop,
nfunc, theta, N, lambda, tau, Tmax, t);
/* 解 u の t=0 におけるグラフを描く */
fmove(a, u[0]);
for (i = 1; i <= N; i++)
fcont(a + i * h, u[i]);
xsync();
/* 係数行列の LU 分解 */
trilu(N - 1, AL + 1, AD + 1, AU + 1);
/* 繰り返し計算 */
nMax = (int) rint(Tmax / tau);
for (n = 1; n <= nMax; n++) {
/* 連立1次方程式を作って解く */
for (i = 1; i < N; i++)
ff[i] = c3 * u[i] + c4 * (u[i - 1] + u[i + 1]);
trisol(N - 1, AL + 1, AD + 1, AU + 1, ff + 1);
for (i = 1; i < N; i++)
u[i] = ff[i];
/* 境界条件 */
u[0] = 0.0;
u[N] = 0.0;
/* 時刻 */
t = n * tau;
if (n % skip == 0) {
/* 数値データの表示 */
if (print_numerical_data)
print100(N, t, u);
/* 解のグラフを描く */
if (erase_always) {
/* これまで描いたものを消し、パラメーターを再表示する */
erase();
print_parameters(xleft, ybottom, xright, ytop,
nfunc, theta, N, lambda, tau, Tmax, t);
}
fmove(a, u[0]);
for (i = 1; i <= N; i++)
fcont(a + i * h, u[i]);
xsync();
}
}
save_picture();
printf("マウスでウィンドウをクリックして下さい。\n");
closepl();
return 0;
}
/**********************************************************/
/* 三重対角行列の LU 分解 (pivoting なし) */
void trilu(int n, double *al, double *ad, double *au)
{
int i, nm1 = n - 1;
/* 前進消去 (forward elimination) */
for (i = 0; i < nm1; i++) {
al[i + 1] /= ad[i];
ad[i + 1] -= au[i] * al[i + 1];
}
}
/* LU 分解済みの三重対角行列を係数に持つ三項方程式を解く */
void trisol(int n, double *al, double *ad, double *au, double *b)
{
int i, nm1 = n - 1;
/* 前進消去 (forward elimination) */
for (i = 0; i < nm1; i++) b[i + 1] -= b[i] * al[i + 1];
/* 後退代入 (backward substitution) */
b[nm1] /= ad[nm1];
for (i = n - 2; i >= 0; i--) b[i] = (b[i] - au[i] * b[i + 1]) / ad[i];
}
/**********************************************************/
/* 初期条件を与える関数。複数登録して番号 nfunc で選択する。 */
double f(double x, int nfunc)
{
/* f(x)=max(x,1-x) */
if (nfunc == 1) {
if (x <= 0.5)
return x;
else
return 1.0 - x;
}
/* f(x)=1 */
else if (nfunc == 2)
return 1.0;
else if (nfunc == 3)
return sin(pi * x);
/* f(x)= 変な形をしたもの(by M.Sakaue) */
else if (nfunc == 4) {
if (x <= 0.1)
return 5 * x;
else if (x <= 0.3)
return -2 * x + 0.7;
else if (x <= 0.5)
return 4.5 * x - 1.25;
else if (x <= 0.7)
return -x + 1.5;
else if (x <= 0.9)
return x + 0.1;
else
return -10 * x + 10.0;
}
/* やはり非対称な形をしたもの(by M.Sakaue) */
else if (nfunc == 5) {
if (x <= 0.2)
return -x + 0.8;
else if (x <= 0.3)
return 4 * x - 0.2;
else if (x <= 0.4)
return -4 * x + 2.2;
else if (x <= 0.7)
return -x + 1.0;
else if (x <= 0.8)
return 4 * x - 2.6;
else
return -x + 1.5;
} else
return 0;
}
/**********************************************************/
/* 配列 u[] の内容(u[0],u[1],...,u[n])を一行 5 個ずつ表示する。 */
void print100(int N, double t, double *u)
{
int i;
printf("T= %12.4e\n", t);
for (i = 0; i < 5; i++)
printf(" I u(i) ");
printf("\n");
for (i = 0; i <= N; i++) {
printf("%3d%12.4e ", i, u[i]);
if (i % 5 == 4)
printf("\n");
}
printf("\n");
}
/**********************************************************/
/* ウィンドウに座標軸を表示する */
void axis(double xleft, double ybottom, double xright, double ytop)
{
linemod("dotted");
fline(xleft, 0.0, xright, 0.0);
fline(0.0, ybottom, 0.0, ytop);
linemod("solid");
}
/* ウィンドウ内の位置 (x,y) から文字列 s を表示する */
void fsymbol(double x, double y, char *s)
{
fmove(x, y);
label(s);
}
#define X(x) (xleft+(x)*(xright-xleft))
#define Y(y) (ybottom+(y)*(ytop-ybottom))
/* パラメーターの値等をウィンドウに表示する */
void print_parameters(double xleft, double ybottom, double xright, double ytop,
int nfunc, double theta, int N,
double lambda, double tau, double Tmax, double t)
{
char message[80];
axis(xleft, ybottom, xright, ytop);
sprintf(message, "heat equation, u(0)=u(1)=0");
fsymbol(X(0.2), Y(0.95), message);
sprintf(message, "nfunc=%d, theta=%g, N=%d, lambda=%g, tau=%g, Tmax=%g",
nfunc, theta, N, lambda, tau, Tmax);
fsymbol(X(0.2), Y(0.9), message);
sprintf(message, "t = %g", t);
if (t != 0.0)
fsymbol(X(0.2), Y(0.85), message);
}
/* 現在ウィンドウに表示されている図をファイルに保存する */
void save_picture(void)
{
char filename[200];
printf("図を保存するファイル名: ");
scanf("%s", filename);
mkplot(filename);
}