32 #include "ope_elementary.h" 38 Tbl _solh_sec_order_r2_pas_prevu (
int,
double,
double,
double,
double,
double,Tbl&) {
40 cout <<
"Homogeneous solution not implemented in Sec_order_r2 : "<< endl ;
52 Tbl _solh_sec_order_r2_r_cheb (
int n,
double alpha,
double beta,
53 double a,
double b,
double c, Tbl& val_lim) {
58 double* coloc =
new double[n] ;
60 int * deg =
new int[3] ;
66 double* air =
new double[n] ;
67 for (
int i=0 ; i<n ; i++)
68 air[i] = alpha*(-
cos(M_PI*i/(n-1))) + beta ;
70 double r_minus = beta-alpha ;
71 double r_plus = beta+alpha ;
74 double delta = (b-a)*(b-a)-4*a*c ;
80 if (fabs(delta) < 1e-14)
87 double lambda_one = ((a-b) +
sqrt(delta)) / 2./a ;
88 double lambda_two = ((a-b) -
sqrt(delta)) / 2./a ;
91 for (
int i=0 ; i<n ; i++)
92 coloc[i] =
pow(air[i], lambda_one) ;
93 cfrcheb(deg, deg, coloc, deg, coloc) ;
94 for (
int i=0 ; i<n ;i++)
95 res.set(0,i) = coloc[i] ;
98 for (
int i=0 ; i<n ; i++)
99 coloc[i] =
pow(air[i], lambda_two) ;
100 cfrcheb(deg, deg, coloc, deg, coloc) ;
101 for (
int i=0 ; i<n ;i++)
102 res.set(1,i) = coloc[i] ;
105 val_lim.set(0,0) =
pow(r_minus, lambda_one) ;
106 val_lim.set(0,1) = lambda_one*
pow(r_minus, lambda_one-1) ;
107 val_lim.set(0,2) =
pow(r_plus, lambda_one) ;
108 val_lim.set(0,3) = lambda_one*
pow(r_plus, lambda_one-1) ;
110 val_lim.set(1,0) =
pow(r_minus, lambda_two) ;
111 val_lim.set(1,1) = lambda_two*
pow(r_minus, lambda_two-1) ;
112 val_lim.set(1,2) =
pow(r_plus, lambda_two) ;
113 val_lim.set(1,3) = lambda_two*
pow(r_plus, lambda_two-1) ;
114 val_lim /=
sqrt(
double(2)) ;
119 double lambda = (a-b)/2./a ;
121 for (
int i=0 ; i<n ; i++)
122 coloc[i] =
pow(air[i], lambda) ;
123 cfrcheb(deg, deg, coloc, deg, coloc) ;
124 for (
int i=0 ; i<n ;i++)
125 res.set(0,i) = coloc[i] ;
128 for (
int i=0 ; i<n ; i++)
129 coloc[i] =
log(air[i])*
pow(air[i], lambda) ;
130 cfrcheb(deg, deg, coloc, deg, coloc) ;
131 for (
int i=0 ; i<n ;i++)
132 res.set(1,i) = coloc[i] ;
135 val_lim.set(0,0) =
pow(r_minus, lambda) ;
136 val_lim.set(0,1) = lambda*
pow(r_minus, lambda-1) ;
137 val_lim.set(0,2) =
pow(r_plus, lambda) ;
138 val_lim.set(0,3) = lambda*
pow(r_plus, lambda-1) ;
140 val_lim.set(1,0) =
log(r_minus)*
pow(r_minus, lambda) ;
141 val_lim.set(1,1) = (1+lambda*
log(r_minus))*
pow(r_minus, lambda-1) ;
142 val_lim.
set(1,2) =
log(r_plus)*
pow(r_plus, lambda) ;
143 val_lim.set(1,3) = (1+lambda*
log(r_plus))*
pow(r_plus, lambda-1) ;
144 val_lim /=
sqrt(
double(2)) ;
149 double real_part = (a-b)/2./a ;
150 double imag_part =
sqrt(-delta)/2./a ;
153 for (
int i=0 ; i<n ; i++)
154 coloc[i] =
pow(air[i], real_part)*
cos(imag_part*
log(air[i])) ;
155 cfrcheb(deg, deg, coloc, deg, coloc) ;
156 for (
int i=0 ; i<n ;i++)
157 res.set(0,i) = coloc[i] ;
160 for (
int i=0 ; i<n ; i++)
161 coloc[i] =
pow(air[i], real_part)*
sin(imag_part*
log(air[i])) ;
162 cfrcheb(deg, deg, coloc, deg, coloc) ;
163 for (
int i=0 ; i<n ;i++)
164 res.set(1,i) = coloc[i] ;
167 val_lim.set(0,0) =
pow(r_minus, real_part)*
cos(imag_part*
log(r_minus)) ;
168 val_lim.set(0,1) = (real_part*
cos(imag_part*
log(r_minus)) -
169 imag_part*
sin(imag_part*
log(r_minus))) *
170 pow(r_minus, real_part-1) ;
171 val_lim.
set(0,2) =
pow(r_plus, real_part)*
cos(imag_part*
log(r_plus)) ;
172 val_lim.set(0,3) = (real_part*
cos(imag_part*
log(r_plus)) -
173 imag_part*
sin(imag_part*
log(r_plus))) *
174 pow(r_plus, real_part-1) ;
176 val_lim.
set(1,0) =
pow(r_minus, real_part)*
sin(imag_part*
log(r_minus)) ;
177 val_lim.set(1,1) = (real_part*
sin(imag_part*
log(r_minus)) +
178 imag_part*
cos(imag_part*
log(r_minus))) *
179 pow(r_minus, real_part-1) ;
180 val_lim.
set(1,2) =
pow(r_plus, real_part)*
sin(imag_part*
log(r_plus)) ;
181 val_lim.set(1,3) = (real_part*
sin(imag_part*
log(r_plus)) +
182 imag_part*
cos(imag_part*
log(r_plus))) *
183 pow(r_plus, real_part-1) ;
184 val_lim /=
sqrt(
double(2)) ;
188 cout <<
"What are you doing here ? Get out or I call the police !" << endl;
204 static Tbl (*solh_sec_order_r2[
MAX_BASE]) (int, double, double,
205 double, double, double,
Tbl&) ;
212 solh_sec_order_r2[i] = _solh_sec_order_r2_pas_prevu ;
215 solh_sec_order_r2[
R_CHEB >>
TRA_R] = _solh_sec_order_r2_r_cheb ;
double alpha
Parameter of the associated mapping.
Cmp log(const Cmp &)
Neperian logarithm.
double s_one_minus
Value of the first homogeneous solution at the inner boundary.
Cmp sqrt(const Cmp &)
Square root.
double ds_two_minus
Value of the derivative of the second homogeneous solution at the inner boundary. ...
double beta
Parameter of the associated mapping.
virtual Tbl get_solh() const
Computes the homogeneous solutions(s).
double ds_two_plus
Value of the derivative of the second homogeneous solution at the outer boundary. ...
double a_param
The parameter a .
Cmp cos(const Cmp &)
Cosine.
int base_r
Radial basis of decomposition.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
double ds_one_plus
Value of the derivative of the first homogeneous solution at the outer boundary.
double b_param
The parameter b .
double s_two_minus
Value of the second homogeneous solution at the inner boundary.
Cmp pow(const Cmp &, int)
Power .
Tbl & set(int l)
Read/write of the value in a given domain.
double s_one_plus
Value of the first homogeneous solution at the outer boundary.
int nr
Number of radial points.
double ds_one_minus
Value of the derivative of the first homogeneous solution at the inner boundary.
Cmp sin(const Cmp &)
Sine.
double c_param
The parameter c .
double s_two_plus
Value of the second homogeneous solution at the outer boundary.
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEB
base de Chebychev ordinaire (fin)