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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
|
#include <iostream>
#include <vector>
#include <cmath>
using namespace std;
using vec = vector< double >;
using matrix = vector< vec >;
//======= Function prototypes ========================================
bool solve( matrix A, vec B, vec &X );
double F( double x, double y, double d ); // The infamous function
//====================================================================
int main()
{
const double epsilon1 = 0.8, epsilon2 = 0.6; // Emissivities
const double sigma = 1.712e-9; // Stefan Boltzmann constant
const double T1 = 1000.0, T2 = 500.0; // Absolute temperatures of the plates
const double d = 1.0; // Separation of the plates
const double w = 1.0; // Width of the plates
for ( int N : { 4, 8, 16, 32, 64, 128, 256 } ) // Number of intervals
{
int NN = 2 * ( N + 1 ); // Number of degrees of freedom (values of u and v)
// Set up coordinates
vec x(1+N);
double dx = w / N;
for ( int i = 0; i <= N; i++ ) x[i] = -0.5 * w + i * dx;
vec y = x;
double dy = dx;
// Set up linear system as (there are other possibilities)
// M {u} = b
// {v}
matrix M( NN, vec(NN,0.0) );
vec B(NN);
vec XX(NN);
double heat1 = epsilon1 * sigma * T1 * T1 * T1 * T1;
double heat2 = epsilon2 * sigma * T2 * T2 * T2 * T2;
for ( int i = 0; i <= N; i++ ) // u equation (first N+1 rows of M; x = x[i])
{
B[i] = heat1;
M[i][i] = 1.0; // u[i] = ...
M[i][N +1] = -( 1.0 - epsilon1 ) * dy * 0.5 * F( x[i], y[0], d ); // "Integral" transferred to LHS
M[i][NN-1] = -( 1.0 - epsilon1 ) * dy * 0.5 * F( x[i], y[N], d );
for ( int j = 1; j < N; j++ ) M[i][N+1+j] = -( 1.0 - epsilon1 ) * dy * F( x[i], y[j], d );
}
for ( int i = 0; i <= N; i++ ) // v equation (last N+1 rows of M; y = y[i])
{
B[N+1+i] = heat2;
M[N+1+i][N+1+i] = 1.0; // v[i] = ....
M[N+1+i][0] = -( 1.0 - epsilon2 ) * dx * 0.5 * F( x[0], y[i], d );
M[N+1+i][N] = -( 1.0 - epsilon2 ) * dx * 0.5 * F( x[N], y[i], d );
for ( int j = 1; j < N; j++ ) M[N+1+i][j] = -( 1.0 - epsilon2 ) * dy * F( x[j], y[i], d );
}
// Solve
cout << "Solving ... " << boolalpha << solve( M, B, XX ) << '\n';
vec U( XX.begin(), XX.begin()+1+N );
vec V( XX.begin()+1+N, XX.end() );
// Post-processing
vec I1(1+N), I2(1+N);
for (int i = 0; i <= N; i++ )
{
I1[i] = ( U[i] - heat1 ) / ( 1.0 - epsilon1 );
I2[i] = ( V[i] - heat2 ) / ( 1.0 - epsilon2 );
}
double Q1 = ( U[0] - I1[0] + U[N] - I1[N] ) * 0.5 * dx;
for ( int i = 1; i < N; i++ ) Q1 += ( U[i] - I1[i] ) * dx;
double Q2 = ( V[0] - I2[0] + V[N] - I2[N] ) * 0.5 * dy;
for ( int i = 1; i < N; i++ ) Q2 += ( V[i] - I2[i] ) * dy;
cout << "N = " << N << " Q1 = " << Q1 << " Q2 = " << Q2 << '\n';
// cout << "\nU values: "; for ( double e : U ) cout << e << " ";
// cout << "\nV values: "; for ( double e : V ) cout << e << " ";
}
}
//====================================================================
double F( double x, double y, double d )
{
double rsq = d * d + ( x - y ) * ( x - y );
return 0.5 / rsq / sqrt( rsq );
}
//====================================================================
bool solve( matrix A, vec B, vec &X )
//--------------------------------------
// Solve AX = B by Gaussian elimination
//--------------------------------------
{
const double SMALL = 1.0e-50;
int n = A.size();
// Row operations for i = 0, ,,,, n - 2 (n-1 not needed)
for ( int i = 0; i < n - 1; i++ )
{
// Pivot: find row r below with largest element in column i
int r = i;
double maxA = abs( A[i][i] );
for ( int k = i + 1; k < n; k++ )
{
double val = abs( A[k][i] );
if ( val > maxA )
{
r = k;
maxA = val;
}
}
if ( r != i )
{
for ( int j = i; j < n; j++ ) swap( A[i][j], A[r][j] );
swap( B[i], B[r] );
}
// Row operations to make upper-triangular
double pivot = A[i][i];
if ( abs( pivot ) < SMALL ) return false; // Singular matrix
for ( int r = i + 1; r < n; r++ ) // On lower rows
{
double multiple = A[r][i] / pivot; // Multiple of row i to clear element in ith column
for ( int j = i; j < n; j++ ) A[r][j] -= multiple * A[i][j];
B[r] -= multiple * B[i];
}
}
if ( abs( A[n-1][n-1] ) < SMALL ) return false; // Singular matrix
// Back-substitute
X = B;
for ( int i = n - 1; i >= 0; i-- )
{
for ( int j = i + 1; j < n; j++ ) X[i] -= A[i][j] * X[j];
X[i] /= A[i][i];
}
return true;
}
//====================================================================
| |
