1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
|
#include <assert.h>
#include <math.h>
#include <iostream>
#include <string.h>
#include <vector>
#include <fstream>
#include <sstream>
#include <stdlib.h>
#include <stdio.h>
/* Subroutines ======================================================
*/
void TA(int N, double *b, double *a, double *c,double *x, double *q)
{
int i;
static double *l=NULL,*u=NULL,*d=NULL,*y=NULL;
if(l == NULL)
{
l = new double[N];
u = new double[N];
d = new double[N];
y = new double[N];
}
/* LU Decomposition */
d[0] = a[0];
u[0] = c[0];
for(i=0;i<N-2;i++)
{
l[i] = b[i]/d[i];
d[i+1] = a[i+1] - l[i]*u[i];
u[i+1] = c[i+1];
}
l[N-2] = b[N-2]/d[N-2];
d[N-1] = a[N-1] - l[N-2]*u[N-2];
/* Forward Substitution [L][y] = [q] */
y[0] = q[0];
for(i=1;i<N;i++)
y[i] = q[i] - l[i-1]*y[i-1];
/* Backward Substitution [U][x] = [y] */
x[N-1] = y[N-1]/d[N-1];
for(i=N-2;i>=0;i--)
x[i] = (y[i] - u[i]*x[i+1])/d[i];
return;
}
void TA_per(int N, double *b, double *a, double *c,double *x, double *q)
{
int i;
double *x1,*x2,*q2;
x1 = new double[N-1];
x2 = new double[N-1];
q2 = new double[N-1];
/* Prepare secondary q */
for(i=0;i<N-1;i++)
q2[i] = 0.0;
q2[0] = -b[N-1];
q2[N-2] = -c[N-2];
TA(N-1,b,a,c,x1,q);
TA(N-1,b,a,c,x2,q2);
x[N-1] = (q[N-1] - c[N-1]*x1[0] - b[N-2]*x1[N-2])/
(a[N-1] + c[N-1]*x2[0] + b[N-2]*x2[N-2]);
for(i=0;i<N-1;i++)
x[i] = x1[i] + x2[i]*x[N-1];
delete[] x1;
delete[] x2;
delete[] q2;
}
void EB(int N, double DN,double *uold, double *unew)
{
TA(N-2,-DN,1.0+2.0*DN,-DN,&unew[1],&uold[1]); ERROR HERE
return;
}
void Cr(int N, double CFL,double *uold, double *unew)
{
double c = 0.25*CFL;
double *q = new double[N];
for(int i=1; i<N-1; i++)
q[i] = uold[i] - 0.25*CFL*(uold[i+1]-uold[i-1]);
q[0] = uold[0] - 0.25*CFL*(uold[1]-uold[N-1]);
q[N-1] = uold[N-1] - 0.25*CFL*(uold[0]-uold[N-2]);
TA_per(N,-c,1.0e0,c,unew,q); AND ERROR HERE
delete[] q;
}
|
