32 #include "ope_elementary.h" 38 Tbl _solh_sec_order_pas_prevu (
int,
double,
double,
double,
double,
double,Tbl&) {
40 cout <<
"Homogeneous solution not implemented in Sec_order : "<< endl ;
52 Tbl _solh_sec_order_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* sigma =
new double[n] ;
67 for (
int i=0 ; i<n ; i++)
68 sigma[i] = alpha*(-
cos(M_PI*i/(n-1))) + beta ;
70 double sigma_minus = beta-alpha ;
71 double sigma_plus = beta+alpha ;
74 double delta = b*b-4*a*c ;
80 if (fabs(delta) < 1e-14)
87 double lambda_one = (-b +
sqrt(delta)) / 2./a ;
88 double lambda_two = (-b -
sqrt(delta)) / 2./a ;
91 for (
int i=0 ; i<n ; i++)
92 coloc[i] =
exp(sigma[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] =
exp(sigma[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) =
exp(sigma_minus*lambda_one) ;
106 val_lim.set(0,1) = lambda_one*
exp(sigma_minus*lambda_one)
108 val_lim.set(0,2) =
exp(sigma_plus*lambda_one) ;
109 val_lim.set(0,3) = lambda_one*
exp(sigma_plus*lambda_one)
112 val_lim.set(1,0) =
exp(sigma_minus*lambda_two) ;
113 val_lim.set(1,1) = lambda_two*
exp(sigma_minus*lambda_two)/
115 val_lim.set(1,2) =
exp(sigma_plus*lambda_two) ;
116 val_lim.set(1,3) = lambda_two*
exp(sigma_plus*lambda_two)
118 val_lim /=
sqrt(
double(2)) ;
123 double lambda = -b/2./a ;
125 for (
int i=0 ; i<n ; i++)
126 coloc[i] =
exp(sigma[i]*lambda) ;
127 cfrcheb(deg, deg, coloc, deg, coloc) ;
128 for (
int i=0 ; i<n ;i++)
129 res.set(0,i) = coloc[i] ;
132 for (
int i=0 ; i<n ; i++)
133 coloc[i] = sigma[i]*
exp(sigma[i]*lambda) ;
134 cfrcheb(deg, deg, coloc, deg, coloc) ;
135 for (
int i=0 ; i<n ;i++)
136 res.set(1,i) = coloc[i] ;
139 val_lim.set(0,0) =
exp(sigma_minus*lambda) ;
140 val_lim.set(0,1) = lambda*
exp(sigma_minus*lambda)/
exp(sigma_minus) ;
141 val_lim.set(0,2) =
exp(sigma_plus*lambda) ;
142 val_lim.set(0,3) = lambda*
exp(sigma_plus*lambda)/
exp(sigma_plus) ;
144 val_lim.set(1,0) = sigma_minus*
exp(sigma_minus*lambda) ;
145 val_lim.set(1,1) =
exp(sigma_minus*lambda)* (lambda*sigma_minus+1)
147 val_lim.
set(1,2) = sigma_plus*
exp(sigma_plus*lambda) ;
148 val_lim.set(1,3) =
exp(sigma_plus*lambda)* (lambda*sigma_plus+1)/
150 val_lim /=
sqrt(
double(2)) ;
155 double real_part = -b/2./a ;
156 double imag_part =
sqrt(-delta)/2./a ;
159 for (
int i=0 ; i<n ; i++)
160 coloc[i] =
exp(sigma[i]*real_part)*
cos(imag_part*sigma[i]) ;
161 cfrcheb(deg, deg, coloc, deg, coloc) ;
162 for (
int i=0 ; i<n ;i++)
163 res.set(0,i) = coloc[i] ;
166 for (
int i=0 ; i<n ; i++)
167 coloc[i] =
exp(sigma[i]*real_part)*
sin(imag_part*sigma[i]) ;
168 cfrcheb(deg, deg, coloc, deg, coloc) ;
169 for (
int i=0 ; i<n ;i++)
170 res.set(1,i) = coloc[i] ;
173 val_lim.set(0,0) =
exp(sigma_minus*real_part)*
cos(imag_part*sigma_minus) ;
174 val_lim.set(0,1) = (real_part*
cos(imag_part*sigma_minus) -
175 imag_part*
sin(imag_part*sigma_minus)) *
exp(real_part*sigma_minus)
177 val_lim.set(0,2) =
exp(sigma_plus*real_part)*
cos(imag_part*sigma_plus) ;
178 val_lim.set(0,3) = (real_part*
cos(imag_part*sigma_plus) -
179 imag_part*
sin(imag_part*sigma_plus)) *
exp(real_part*sigma_plus)
182 val_lim.set(1,0) =
exp(sigma_minus*real_part)*
sin(imag_part*sigma_minus) ;
183 val_lim.set(1,1) = (real_part*
sin(imag_part*sigma_minus) +
184 imag_part*
cos(imag_part*sigma_minus)) *
exp(real_part*sigma_minus)
186 val_lim.set(1,2) =
exp(sigma_plus*real_part)*
sin(imag_part*sigma_plus);
187 val_lim.set(1,3) = (real_part*
sin(imag_part*sigma_plus) +
188 imag_part*
cos(imag_part*sigma_plus)) *
exp(real_part*sigma_plus)
190 val_lim /=
sqrt(
double(2)) ;
194 cout <<
"What are you doing here ? Get out or I call the police !" << endl;
210 static Tbl (*solh_sec_order[
MAX_BASE]) (int, double, double,
211 double, double, double,
Tbl&) ;
218 solh_sec_order[i] = _solh_sec_order_pas_prevu ;
221 solh_sec_order[
R_CHEB >>
TRA_R] = _solh_sec_order_r_cheb ;
double alpha
Parameter of the associated mapping.
Cmp exp(const Cmp &)
Exponential.
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.
double ds_two_plus
Value of the derivative of the second homogeneous solution at the outer boundary. ...
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 b_param
The parameter .
double ds_one_plus
Value of the derivative of the first homogeneous solution at the outer boundary.
double c_param
The parameter .
double a_param
The parameter .
double s_two_minus
Value of the second homogeneous solution at the inner boundary.
virtual Tbl get_solh() const
Computes the homogeneous solutions(s).
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 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)