авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 | 2 ||

«ПЕРМСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ УДК 517.97+517.688 На правах рукописи ...»

-- [ Страница 3 ] --

} DOUBLE& h i ( ) { return u ;

} INTERVAL& ( const INTERVAL& x);

operator+= INTERVAL& ( const INTERVAL& x);

operator= INTERVAL& ( const INTERVAL& x);

operator= INTERVAL& ( const INTERVAL& x);

operator/= };

MATH: : INTERVAL SQRT ( const MATH: : INTERVAL& x ) ;

};

i n l i n e MATH: : INTERVAL operator+ ( const MATH: : INTERVAL& x, MATH: : INTERVAL& y ) const { MATH: : INTERVAL z = x ;

z += y ;

return z ;

} i n l i n e MATH: : INTERVAL operator ( const MATH: : INTERVAL& x, MATH: : INTERVAL& y ) const { MATH: : INTERVAL z = x ;

z = y ;

return z ;

} i n l i n e MATH: : INTERVAL operator ( const MATH: : INTERVAL& x, MATH: : INTERVAL& y ) const { MATH: : INTERVAL z = x ;

z = y ;

return z ;

} i n l i n e MATH: : INTERVAL operator / ( const MATH: : INTERVAL& x, MATH: : INTERVAL& y ) const { MATH: : INTERVAL z = x ;

z /= y ;

return z ;

} // При выполнении операций сравнения предполагается, // что интервалы нельзя назвать "равными"только в том случае, // если они не пересекаются [a1, a2 ] [b1, b2 ] = i n l i n e bool operator ( const MATH: : INTERVAL& x, const MATH: : INTERVAL& y ) { return x. h i ( ) y. l o ( ) ;

} i n l i n e bool operator ( const MATH: : INTERVAL& x, const MATH: : INTERVAL& y ) { return x. l o ( ) y. h i ( ) ;

} // Более мягкая операция сравнения x = [x1, x2 ] и y = [y1, y2 ] // -2 : x2 y // -1 : x1 y1 и x2 y // 0 : x y или y x // 1 : x1 y1 и x2 y // 2 : x1 y i n t compare ( const MATH: : INTERVAL& x, const MATH: : INTERVAL& y ) ;

MATH: : INTERVAL ABS ( const MATH: : INTERVAL& x ) ;

s t d : : ostream& operator ( s t d : : ostream& os, const MATH: : INTERVAL& x ) ;

#e ndif interval.cpp #i f d e f _WIN #include f l o a t. h #e l s e #include f e n v. h #e ndif #include mpfr. h #include cmath #include " i n t e r v a l. hpp" namespace MATH { #i f d e f _WIN typedef unsigned i n t FP_STATUS;

#e l s e typedef i n t FP_STATUS ;

#e ndif // режим округления вниз i n l i n e FP_STATUS round_down ( ) { #i f d e f _WIN return c o n t r o l _ f p (_MCW_RC,_RC_DOWN) ;

#e l s e return f e s e t r o u n d (FE_DOWNWARD) ;

#e ndif } // режим округления вверх i n l i n e FP_STATUS round_up ( ) { #i f d e f _WIN return c o n t r o l _ f p (_MCW_RC,_RC_UP) ;

#e l s e return f e s e t r o u n d (FE_UPWARD) ;

#e ndif } // задание/восстановление направления округления i n l i n e FP_STATUS s e t _ r c _ s t a t u s (FP_STATUS s t a t u s ) { #i f d e f _WIN return c o n t r o l _ f p (_MCW_RC s t a t u s ) ;

, #e l s e return f e s e t r o u n d ( s t a t u s ) ;

#e ndif } INTERVAL : : INTERVAL ( mpq_class& x ) { mpfr_t x_mp;

m p f r _ i n i t 2 (x_mp, s i z e o f (MATH: : DOUBLE) 1);

mpfr_set_q (x_mp, x. get_mpq_t ( ),GMP_RNDD) ;

l = mpfr_get_ld (x_mp,GMP_RNDD) ;

m p f r _ c l e a r (x_mp ) ;

m p f r _ i n i t 2 (x_mp, s i z e o f (MATH: : DOUBLE) 1);

mpfr_set_q (x_mp, x. get_mpq_t ( ),GMP_RNDU) ;

u = mpfr_get_ld (x_mp,GMP_RNDU) ;

m p f r _ c l e a r (x_mp ) ;

} INTERVAL : : INTERVAL ( long numer, long denom ) { FP_STATUS prev_mode = round_down ( ) ;

l = (DOUBLE) numer / (DOUBLE) denom ;

round_up ( ) ;

u = (DOUBLE) numer / (DOUBLE) denom ;

s e t _ r c _ s t a t u s ( prev_mode ) ;

} INTERVAL& INTERVAL : : operator+= ( const MATH: : INTERVAL& x ) { FP_STATUS prev_mode = round_down ( ) ;

l += x. l ;

round_up ( ) ;

u += x. u ;

s e t _ r c _ s t a t u s ( prev_mode ) ;

return t h i s ;

} INTERVAL& INTERVAL : : operator= ( const MATH: : INTERVAL& x ) { FP_STATUS prev_mode = round_down ( ) ;

l = x. u ;

round_up ( ) ;

u = x. l ;

s e t _ r c _ s t a t u s ( prev_mode ) ;

return t h i s ;

} INTERVAL& INTERVAL : : operator= ( const MATH: : INTERVAL& x ) { INTERVAL r e s ;

DOUBLE z [ 4 ] ;

FP_STATUS prev_mode = round_down ( ) ;

z [0] = l x. l ;

z [1] = l x.u;

z [2] = u x. l ;

z [3] = u x.u;

DOUBLE L = MIN(MIN( z [ 0 ], z [ 1 ] ), MIN( z [ 2 ], z [ 3 ] ) ) ;

round_up ( ) ;

z [0] = l x. l ;

z [1] = l x.u;

z [2] = u x. l ;

z [3] = u x.u;

u = MAX(MAX( z [ 0 ], z [ 1 ] ),MAX( z [ 2 ], z [ 3 ] ) ) ;

l = L;

s e t _ r c _ s t a t u s ( prev_mode ) ;

return t h i s ;

} INTERVAL& INTERVAL : : operator/= ( const MATH: : INTERVAL& x ) { i f ( x. l =0 && x. u=0) throw "INTERVAL ( / ) : делениенаноль" ;

INTERVAL r e s ;

DOUBLE z [ 4 ] ;

FP_STATUS prev_mode = round_down ( ) ;

z [0] = l / x. l ;

z [1] = l / x.u;

z [2] = u / x. l ;

z [3] = u / x.u;

DOUBLE L = MIN(MIN( z [ 0 ], z [ 1 ] ), MIN( z [ 2 ], z [ 3 ] ) ) ;

round_up ( ) ;

z [0] = l / x. l ;

z [1] = l / x.u;

z [2] = u / x. l ;

z [3] = u / x.u;

u = MAX(MAX( z [ 0 ], z [ 1 ] ),MAX( z [ 2 ], z [ 3 ] ) ) ;

l= L ;

s e t _ r c _ s t a t u s ( prev_mode ) ;

return t h i s ;

} MATH: : INTERVAL SQRT ( const MATH: : INTERVAL& x ) { MATH: : INTERVAL r e s ;

FP_STATUS prev_mode = round_down ( ) ;

res. lo () = sqrt (x. lo ( ) ) ;

round_up ( ) ;

res. hi () = sqrt (x. hi ( ) ) ;

s e t _ r c _ s t a t u s ( prev_mode ) ;

return r e s ;

} };

MATH: : INTERVAL ABS ( const MATH: : INTERVAL& x ) { i f ( x. l o ()=0) return x ;

else i f ( x. h i ()=0) return MATH: : INTERVAL(x. h i (), x. l o ( ) ) ;

e l s e return MATH: : INTERVAL( 0,MATH: :MAX(x. l o ( ), x. h i ( ) ) ) ;

} i n t compare ( const MATH: : INTERVAL& x, const MATH: : INTERVAL& y ) { i f ( x. h i ( ) y. l o ( ) ) return 2;

e l s e i f ( ( x. l o () y. l o ( ) ) && ( x. h i () y. h i ( ) ) ) return 1;

e l s e i f ( x. l o ( ) y. h i ( ) ) return 2 ;

e l s e i f ( ( x. l o () y. l o ( ) ) && ( x. h i () y. h i ( ) ) ) return 1 ;

e l s e return 0 ;

} s t d : : ostream& operator ( s t d : : ostream& os, const MATH: : INTERVAL& x ) { i f ( x. l o ()==x. h i ( ) ) o s x. l o ( ) ;

else o s " [ " x. l o ( ) " ;

" x. h i ( ) " ] " ;

return o s ;

} A.4 Постановка задачи problem.hpp #i f n d e f PROBLEM #define PROBLEM v e c t o r #include v a l a r r a y #include s t r i n g #include gmpxx. h #include " l i n a l g. hpp" #include " d a t a u t i l. hpp" #include " o r t h o. hpp" #include " i n t e r v a l. hpp" #include #define _PROB_ #define PROBLEM PROBLEM_ c l a s s F_LAMBDA : public MATH: : PARAMETERS { public :

s i z e _ t number ;

F_LAMBDA ( s i z e _ t num = 0 ) : number (num) {} };

namespace PROBLEM { // границы области extern mpq_class a, b ;

// число подынтервалов extern s i z e _ t n ;

// границы подынтервалов extern s t d : : v a l a r r a y mpq_class tau ;

// подрядок аппроксимации на подынтервале extern s t d : : v a l a r r a y s i z e _ t m;

// правая часть и ядро интегрального уравнения MATH: : DOUBLE f ( const s t d : : v a l a r r a y MATH: : DOUBLE& t, MATH: : PARAMETERS p = NULL ) ;

MATH: : DOUBLE K ( const s t d : : v a l a r r a y MATH: : DOUBLE& t, MATH: : PARAMETERS p = NULL ) ;

// используемое множество ортогональных функций typedef MATH: : LEGENDRE BASIS_FN ;

// рациональные typedef MATH: : LEGENDRE_DBL BASIS_FN_DBL;

// DOUBLE-представление // название задачи (для определения имён файлов) extern s t d : : s t r i n g name ;

// инициализация и очистка void i n i t ( ) ;

void d e s t r o y ( ) ;

//================================================ // множители Лагранжа extern s t d : : v a l a r r a y mpq_class lambda ;

//================================================ // область и моном коэффициента ряда Фурье функции st ruct COEFF_RD { // индексы подобласти и мультистепени size_t region, degree ;

};

extern s t d : : v e c t o r COEFF_RD c o e f f _ r d ;

// значения коэффициентов для f extern s t d : : v e c t o r s t d : : v a l a r r a y mpq_class f _ c o e f f ;

// значения коэффициентов для K extern s t d : : v e c t o r s t d : : v a l a r r a y mpq_class K_coeff ;

// базисные функции extern s t d : : v e c t o r MATH: : ORTHOmpq_class b a s i s _ f n ;

extern s t d : : v e c t o r MATH: : ORTHO MATH: : DOUBLE b a s i s _ f n _ d b l ;

// вычисление коэффициентов рядов Фурье для f и K void c r e a t e _ c o e f f s ( ) ;

// вычисление значений аппроксимирующих функций MATH: : DOUBLE K_approx ( s t d : : v a l a r r a y MATH: : DOUBLE& t ) ;

MATH: : DOUBLE f_approx ( s t d : : v a l a r r a y MATH: : DOUBLE& t ) ;

//================================================ // СЛАУ extern MATH: : MATRIXmpq_class A, B, Z, invA ;

// создание матрицы коэффициентов СЛАУ и матрицы (вектора) правой части void create_matrix_A ( ) ;

void create_ v ector_ B ( ) ;

// решение СЛАУ и вычисление матрицы, обратной к A // вычисление коэффициентов решения интегрального уравнения (Z) void solve_IEF ( ) ;

// вычисление значений множителя Лагранжа void s o l v e _ L a g r a n g e ( ) ;

// оценка квадрата нормы резольвентного оператора MATH: : INTERVAL estimate_N2RC ( ) ;

// оценка верхней границы спектра K void ub_spectrum ( mpq_class X0 ) ;

// кусочно–полиномиальное представление решения // уравнения Фредгольма (z) и вариационной задачи (x) extern s t d : : v e c t o r MATH: : POLYNOM MATH: : DOUBLE z, x ;

//================================================ // Значение решения интегрального уравнения в точке t MATH: : DOUBLE Z s o l u t i o n ( s t d : : v a l a r r a y MATH: : DOUBLE& t ) ;

//================================================ // Значение решения вариационной задачи в точке t MATH: : DOUBLE X s o l u t i o n ( s t d : : v a l a r r a y MATH: : DOUBLE& t ) ;

//================================================ // z extern mpq_class normZ2, normK2 ;

};

#e ndif problem.cpp // для использования длинного представления // математических констант // ( определяется как M_PI) #define USE_GNU #include cmath #include c s t d i o #include " i n t e g r a l. hpp" #include " problem. hpp" namespace PROBLEM { //======================================== // данные, а также функции инициализации // и очистки, одинаковые для всех задач s t d : : v e c t o r COEFF_RD c o e f f _ r d ;

s t d : : v e c t o r s t d : : v a l a r r a y mpq_class f _ c o e f f ;

s t d : : v e c t o r s t d : : v a l a r r a y mpq_class K_coeff ;

std : : vectorMATH: : ORTHOmpq_class b a s i s _ f n ;

std : : vectorMATH: : ORTHO MATH: : DOUBLE b a s i s _ f n _ d b l ;

std : : vectorMATH: : POLYNOM MATH: : DOUBLE z, x ;

s i z e _ t N,M;

void main_ init ( ) { // число коэффициентов (размерность СЛАУ) N = m[ 0 ] + 1 ;

f o r ( s i z e _ t i = 1 ;

i n ;

i ++) N += m[ i ] + 1 ;

// определяем максимальную степень аппроксимации s i z e _ t M = m. max ( ) + 1 ;

// определяем комбинации индексов областей и степеней c o e f f _ r d. r e s i z e (N ) ;

size_t idx = 0;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { f o r ( s i z e _ t j = 0 ;

j = m[ i ] ;

j ++) { coeff_rd [ idx ]. region = i ;

coeff_rd [ idx ]. degree = j ;

i d x ++;

} } // генерируем базисные функции b a s i s _ f n. r e s i z e (M) ;

b a s i s _ f n _ d b l. r e s i z e (M) ;

f o r ( s i z e _ t i = 0 ;

i b a s i s _ f n. s i z e ( ) ;

i ++) { b a s i s _ f n [ i ] = new BASIS_FN( i ) ;

b a s i s _ f n _ d b l [ i ] = new BASIS_FN_DBL( i ) ;

} } void main_destroy ( ) { // удаляем базисные функции for ( size_t i = 0 ;

i basis_fn. s i z e ( ) ;

i ++) { de le t e b a s i s _ f n [ i ] ;

de le t e b a s i s _ f n _ d b l [ i ] ;

} } };

//======================================== namespace PROBLEM_ { mpq_class a = mpq_class ( 0 ), b = mpq_class ( 1 0 ) ;

mpq_class a l p h a 1 = mpq_class ( 1);

// mpq_class a l p h a 2 = mpq_class ( 1 ) ;

// size_t n ;

s t d : : v a l a r r a y s i z e _ t m;

s t d : : v a l a r r a y mpq_class tau ;

s t d : : s t r i n g name ( " data /ex_1" ) ;

s t d : : v a l a r r a y mpq_class lambda ( mpq_class ( 0 ), 1 ) ;

// Первообразная функции p(t) t // P (t) = 101 e 5 10 cos(2t) + sin(2t) MATH: : DOUBLE P (MATH: : DOUBLE t ) { using namespace MATH;

s t a t i c DOUBLE c = DOUBLE( 5)/DOUBLE( 1 0 1 ) ;

DOUBLE t 2 = DOUBLE( 2 ) t ;

return c exp( t /DOUBLE( 5 ) ) (DOUBLE( 1 0 ) c o s ( t 2 )+ s i n ( t 2 ) ) ;

} // Первообразная функции q(t) // Q(t) = t + sin(t) MATH: : DOUBLE Q (MATH: : DOUBLE t ) { return t+s i n ( t ) ;

} // Правая часть уравнения Фредгольма f + f · Q(t)Q(b) // f = 1 P (t) P (b) + MATH: : DOUBLE f ( const s t d : : v a l a r r a y MATH: : DOUBLE& t, MATH: : PARAMETERS p ) { i f ( p == NULL) throw "Неопределённомер lambda \n" ;

MATH: : DOUBLE r e s = 0 ;

s t a t i c MATH: : DOUBLE B = b. get_d ( ) ;

s t a t i c MATH: : DOUBLE ALPHA1 = a l p h a 1. get_d ( ) ;

switch ( static_castF_LAMBDA(p)number ) { case 0 :

r e s = ALPHA1 (P( t [ 0 ] ) P(B) ) + (Q( t [ 0 ] ) Q(B) ) /MATH: : DOUBLE( 2 ) ;

break ;

case 1 :

r e s = 0.5;

break ;

de fault :

throw "Неправильныйномер lambda \n" ;

} return r e s ;

} // Ядро уравнения Фредгольма // K(t, s) = P max(t, s) P (b) MATH: : DOUBLE K ( const s t d : : v a l a r r a y MATH: : DOUBLE& t, MATH: : PARAMETERS ) { s t a t i c MATH: : DOUBLE B = b. get_d ( ) ;

return P(MATH: :MAX( t [ 0 ], t [ 1 ] ) ) P(B ) ;

} void i n i t ( ) { n = 15;

m. r e s i z e ( n ) ;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) m[ i ] = 3 ;

char b u f f [ 4 0 ] ;

s p r i n t f ( b u f f, "_%lu_%l u ", n,m [ 0 ] ) ;

name += b u f f ;

tau. r e s i z e ( n + 1);

mpq_class h = ( ba ) / mpq_class ( n ) ;

tau [ 0 ] = a ;

f o r ( s i z e _ t i = 1 ;

i = n ;

i ++) tau [ i ] = tau [ i 1]+h ;

#i f d e f _PROB_ main_ init ( ) ;

#e ndif } // init void d e s t r o y ( ) { #i f d e f _PROB_ main_destroy ( ) ;

#e ndif } };

//======================================== namespace PROBLEM_ { mpq_class a = 0, b = 1 0 ;

mpq_class a l p h a 1 = 1, a l p h a 2 = 1 ;

mpq_class b e t a 1 = 2, b e t a 2 = 3 ;

size_t n ;

s t d : : v a l a r r a y s i z e _ t m;

s t d : : v a l a r r a y mpq_class tau ;

s t d : : s t r i n g name ( " data /ex_2" ) ;

s t d : : v a l a r r a y mpq_class lambda ( mpq_class ( 0 ), 1 ) ;

// начальная функция 1 (t) MATH: : DOUBLE phi_1 (MATH: : DOUBLE t ) { return 1 ;

} // конечная функция 2 (t) MATH: : DOUBLE phi_2 (MATH: : DOUBLE t ) { return 1;

} // Первообразная для p(t) t // P (t) = 101 e 5 10 cos(2t) + sin(2t) MATH: : DOUBLE P (MATH: : DOUBLE t ) { using namespace MATH;

s t a t i c MATH: : DOUBLE c = DOUBLE( 5)/DOUBLE( 1 0 1 ) exp( t /DOUBLE( 5 ) ) ;

return c (DOUBLE( 1 0 ) c o s (DOUBLE( 2 ) t )+ s i n (DOUBLE( 2 ) t ) ) ;

} // Первообразная для q(t) // Q(t) = sin(t) + t MATH: : DOUBLE Q (MATH: : DOUBLE t ) { return s i n ( t )+ t ;

} // Первообразная для g(t) // G(t) = 5 arctg(t) MATH: : DOUBLE G (MATH: : DOUBLE t ) { return atan ( t ) MATH: : DOUBLE( 5 ) ;

} // g(t)1 (t 1 ) MATH: : DOUBLE GF1 (MATH: : DOUBLE t ) { return atan ( t ) MATH: : DOUBLE( 5 ) ;

} // g(t)2 (t + 2 ) MATH: : DOUBLE GF2 (MATH: : DOUBLE t ) { return atan ( t ) MATH: : DOUBLE( 5 ) ;

} // характеристическая функция для [a, b] MATH: : DOUBLE c h i (MATH: : DOUBLE t, MATH: : DOUBLE a, MATH: : DOUBLE b ) { return ( a=t && t=b ) ? MATH: : DOUBLE( 1 ) : 0 ;

} // функции отклонения аргумента MATH: : DOUBLE h1 (MATH: : DOUBLE t) { s t a t i c MATH: : DOUBLE BETA1 = b e t a 1. get_d ( ) ;

return tBETA1;

} MATH: : DOUBLE h2 (MATH: : DOUBLE t) { s t a t i c MATH: : DOUBLE BETA2 = b e t a 2. get_d ( ) ;

return t+BETA2;

} // правая часть MATH: : DOUBLE f ( const s t d : : v a l a r r a y MATH: : DOUBLE& t, MATH: : PARAMETERS p ) { i f ( p == NULL) throw "Неопределённомер lambda \n" ;

using namespace MATH;

s t a t i c MATH: : DOUBLE ALPHA1 = a l p h a 1. get_d ( ) ;

s t a t i c MATH: : DOUBLE BETA1 = b e t a 1. get_d ( ), BETA2 = b e t a 2. get_d ( ) ;

s t a t i c MATH: : DOUBLE A = a. get_d ( ), B = b. get_d ( ) ;

MATH: : DOUBLE r e s = 0 ;

switch ( static_cast F_LAMBDA(p)number ) { case 0 :

{ DOUBLE f 0 = DOUBLE( 2 ) ALPHA1 (P(B)P( t [ 0 ] ) ) + Q(B)Q( t [ 0 ] ) ;

DOUBLE f 1 = GF1(MIN(A+BETA1, BBETA2) ) GF1(MIN(MAX(A, t [0] BETA2),MIN(A+BETA1, BBETA2) ) ) + ALPHA1 (G(BBETA2)G(MAX(A+BETA1, t [0] BETA2 ) ) ) ;

DOUBLE f 2 = ALPHA1 (G(MAX(BBETA2, t [ 0 ] +BETA1) ) G( t [ 0 ] +BETA1) ) + GF2(MAX(B, t [ 0 ] +BETA1) ) GF2(MAX(BBETA2, t [ 0 ] +BETA1 ) ) ;

r e s = ( f 0+f 1+f 2 ) /DOUBLE( 2 ) ;

} break ;

case 1 :

r e s = 0.5;

break ;

de fault :

throw "Неправильныйномер lambda \n" ;

} return r e s ;

} // ядро MATH: : DOUBLE K ( const s t d : : v a l a r r a y MATH: : DOUBLE& t, MATH: : PARAMETERS ) { using namespace MATH;

s t a t i c MATH: : DOUBLE BETA1 = b e t a 1. get_d ( ), BETA2 = b e t a 2. get_d ( ) ;

s t a t i c MATH: : DOUBLE B = b. get_d ( ) ;

DOUBLE K0 = P(B) P(MAX( t [ 0 ], t [ 1 ] ) ) ;

DOUBLE a1 = MIN(BBETA2,MAX( t [ 0 ] +BETA1, t [1] BETA2 ) ) ;

DOUBLE a2 = MIN(BBETA2,MAX( t [ 1 ] +BETA1, t [0] BETA2 ) ) ;

DOUBLE K1 = G(BBETA2) (G( a1)+G( a2 ) ) /DOUBLE( 2 ) ;

return (K0+K1 ) ;

} void i n i t ( ) { n = 10;

m. r e s i z e ( n ) ;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) m[ i ] = 3 ;

char b u f f [ 4 0 ] ;

s p r i n t f ( b u f f, "_%lu_%l u ", n,m [ 0 ] ) ;

name += b u f f ;

tau. r e s i z e ( n + 1);

mpq_class h = ( ba ) / mpq_class ( n ) ;

tau [ 0 ] = a ;

f o r ( s i z e _ t i = 1 ;

i = n ;

i ++) tau [ i ] = tau [ i 1] + h ;

#i f d e f _PROB_ main_ init ( ) ;

#e ndif } // init void d e s t r o y ( ) { #i f d e f _PROB_ main_destroy ( ) ;

#e ndif } };

A.5 Доказательный вычислительный экспе римент fredholm.cpp #include i o s t r e a m #include f s t r e a m #include " problem. hpp" #include " approx. hpp" #include " i n t e g r a l. hpp" namespace PROBLEM { void c r e a t e _ c o e f f s ( ) { // параметры для аппроксимации std : : valarray MATH: : DOUBLE l b ( 1 ), ub ( 1 ) ;

std : : vector MATH: : ORTHO MATH: : DOUBLE o r t h o ( 1 ) ;

// аппроксимация f s t d : : c o u t "Аппроксимация f... " ;

std : : cout. f l u s h ( ) ;

f _ c o e f f. r e s i z e ( lambda. s i z e ( ) + 1 ) ;

MATH: : PARAMETERS LAMBDA = new F_LAMBDA;

f o r ( s i z e _ t L = 0 ;

L = lambda. s i z e ( ) ;

L++) { f_coeff [L ]. r e s i z e ( coeff_rd. s i z e ( ) ) ;

f o r ( s i z e _ t i = 0 ;

i f _ c o e f f [ L ]. s i z e ( ) ;

i ++) { // границы области аппроксимации l b [ 0 ] = tau [ c o e f f _ r d [ i ]. r e g i o n ]. get_d ( ) ;

ub [ 0 ] = tau [ c o e f f _ r d [ i ]. r e g i o n + 1 ]. get_d ( ) ;

// в ortho содержатся указатели на базисные функции, // по которым вычисляется коэффициент ряда Фурье // degreek ortho [ 0 ] = basis_fn_dbl [ coeff_ rd [ i ]. degree ] ;

// вычисляем коэффициент F_LAMBDA(LAMBDA)number = L ;

static_cast f_coeff [L ] [ i ] = mpq_class ( double (MATH: : f o u r i e r _ c o e f f (PROBLEM: : f, lb, ub, ortho, 1 e 5,LAMBDA) ) ) ;

} } de le t e LAMBDA;

s t d : : c o u t "Ok\n" ;

// аппроксимация K s t d : : c o u t "АппроксимацияK... " ;

std : : cout. f l u s h ( ) ;

lb. r e s i z e ( 2 ) ;

ub. r e s i z e ( 2 ) ;

ortho. r e s i z e ( 2 ) ;

K_coeff. r e s i z e ( c o e f f _ r d. s i z e ( ) ) ;

f o r ( s i z e _ t i = 0 ;

i K_coeff. s i z e ( ) ;

i ++) { K_coeff [ i ]. r e s i z e ( c o e f f _ r d. s i z e ( ) ) ;

// границы области аппроксимации // (первый набор) l b [ 0 ] = tau [ c o e f f _ r d [ i ]. r e g i o n ]. get_d ( ) ;

ub [ 0 ] = tau [ c o e f f _ r d [ i ]. r e g i o n + 1 ]. get_d ( ) ;

// базисные функции (первый набор) ortho [ 0 ] = basis_fn_dbl [ coeff_ rd [ i ]. degree ] ;

f o r ( s i z e _ t j = 0 ;

j K_coeff [ i ]. s i z e ( ) ;

j ++) { // границы области аппроксимации // (второй набор) l b [ 1 ] = tau [ c o e f f _ r d [ j ]. r e g i o n ]. get_d ( ) ;

ub [ 1 ] = tau [ c o e f f _ r d [ j ]. r e g i o n + 1 ]. get_d ( ) ;

// базисные функции (второй набор) ortho [ 1 ] = basis_fn_dbl [ coeff_ rd [ j ]. degree ] ;

// вычисляем коэффициент K_coeff [ i ] [ j ] = mpq_class ( double (MATH: : f o u r i e r _ c o e f f (PROBLEM: : K, lb, ub, o r t h o ) ) ) ;

} } s t d : : c o u t "Ok\n" ;

} //================================== MATH: : MATRIX mpq_class A;

void create_matrix_A ( ) { s t d : : c o u t "СозданиематрицыA... " ;

std : : cout. f l u s h ( ) ;

// размерность СЛАУ size_t N = coeff_rd. s i z e ( ) ;

A. r e s i z e (N, N ) ;

// задаём элементы матрицы в зависимости от подобластей // и мультииндексов коэффициентов K f o r ( s i z e _ t i = 0 ;

i N;

i ++) f o r ( s i z e _ t j = 0 ;

j N;

j ++) { A( i, j ) = i==j ? mpq_class ( 1 ) : mpq_class ( 0 ) ;

mpq_class v a l = K_coeff [ i ] [ j ] ;

v a l = ( tau [ c o e f f _ r d [ j ]. r e g i o n +1] tau [ c o e f f _ r d [ j ]. r e g i o n ] ) ;

mpq_class o r t h o = b a s i s _ f n [ c o e f f _ r d [ j ]. d e g r e e ] ;

MATH: : ORTHO v a l = ortho norm2 ( ) / ( ortho h i ( ) ortho l o ( ) ) ;

A( i, j ) = v a l ;

} s t d : : c o u t "Ok\n" ;

} MATH: : MATRIXmpq_class B ;

void create_ v ector_ B ( ) { s t d : : c o u t "СозданиематрицыB... " ;

std : : cout. f l u s h ( ) ;

// размерность СЛАУ size_t N = coeff_rd. s i z e ( ) ;

B. r e s i z e (N, lambda. s i z e ( ) + 1 ) ;

// задаём элементы матрицы в зависимости от подобластей // и мультииндексов коэффициентов K f o r ( s i z e _ t i = 0 ;

i N;

i ++) { f o r ( s i z e _ t L = 0 ;

L = lambda. s i z e ( ) ;

L++) { B( i, L) = mpq_class ( 0 ) ;

} f o r ( s i z e _ t j = 0 ;

j N;

j ++) { mpq_class d = ( tau [ c o e f f _ r d [ j ]. r e g i o n +1] tau [ c o e f f _ r d [ j ]. r e g i o n ] ) ;

mpq_class o r t h o = b a s i s _ f n [ c o e f f _ r d [ j ]. d e g r e e ] ;

MATH: : ORTHO d = ortho norm2 ( ) / ( ortho h i ( ) ortho l o ( ) ) ;

f o r ( s i z e _ t L = 0 ;

L = lambda. s i z e ( ) ;

L++) { B( i, L) += d K_coeff [ i ] [ j ] f _ c o e f f [ L ] [ j ] ;

} } } s t d : : c o u t "Ok\n" ;

} //================================== MATH: : MATRIX MATH: : POLYNOM mpq_class z_polynom_lambda ;

MATH: : MATRIX mpq_class invA, Z ;

void solve_IEF ( ) { s t d : : c o u t "ОбращениематрицыA... " ;

std : : cout. f l u s h ( ) ;

invA = i n v e r s e (A ) ;

s t d : : c o u t "Ok\РешениеnСЛАУ... " ;

std : : cout. f l u s h ( ) ;

Z = invA B ;

// кусочно-полиномиальное представление решения s t d : : c o u t "Ok" s t d : : e n d l "Кусочнополиномиальноепредставлениерешения... " ;

std : : cout. f l u s h ( ) ;

z_polynom_lambda. r e s i z e ( n, lambda. s i z e ( ) + 1 ) ;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) f o r ( s i z e _ t L = 0 ;

L = lambda. s i z e ( ) ;

L++) { z_polynom_lambda ( i, L ). r e s i z e ( 1 ) ;

z_polynom_lambda ( i, L ) [ 0 ] = 0 ;

} // z(t) = N ui (t) L l (fil + Zil ) i=1 l= // или // z(t) = L l N ui (t)(fil + Zil ) l=0 i= f o r ( s i z e _ t i = 0 ;

i Z. rows ( ) ;

i ++) f o r ( s i z e _ t L = 0 ;

L = lambda. s i z e ( ) ;

L++) { mpq_class c o e f f = f _ c o e f f [ L ] [ i ] + Z ( i, L ) ;

z_polynom_lambda ( c o e f f _ r d [ i ]. r e g i o n, L ) += c o e f f ( dynamic_cast MATH: : POLYNOM mpq_class ( basis_fn [ coeff_rd [ i ]. degree ] ) ) ;

} s t d : : c o u t "Ok\n" ;

} //================================== MATH: : INTERVAL estimate_N2RC ( ) { // размерность СЛАУ size_t N = coeff_rd. s i z e ( ) ;

// предварительные преобразования к интервальному виду MATH: : MATRIX MATH: : INTERVAL invA_int ;

invA_int. r e s i z e (N,N ) ;

f o r ( s i z e _ t i = 0 ;

i N;

i ++) f o r ( s i z e _ t j = 0 ;

j N;

j ++) { invA_int ( i, j ) = MATH: : INTERVAL( invA ( i, j ) ) ;

} std : : valarray MATH: : INTERVAL tau_ int ( tau. s i z e ( ) ) ;

f o r ( s i z e _ t i = 0 ;

i tau. s i z e ( ) ;

i ++) { tau_ int [ i ] = MATH: : INTERVAL( tau [ i ] ) ;

} std : : valarray MATH: : INTERVAL norm2_int ( b a s i s _ f n. s i z e ( ) ) ;

f o r ( s i z e _ t i = 0 ;

i b a s i s _ f n. s i z e ( ) ;

i ++) { mpq_class nrm = b a s i s _ f n [ i ]norm2 ( ) ;

norm2_int [ i ] = MATH: : INTERVAL( nrm ) ;

} MATH: : MATRIX MATH: : INTERVAL K _ coeff_ int ;

K _ coeff_ int. r e s i z e (N,N ) ;

f o r ( s i z e _ t i = 0 ;

i N;

i ++) f o r ( s i z e _ t j = 0 ;

j N;

j ++) { K _ coeff_ int ( i, j ) = MATH: : INTERVAL( K_coeff [ i ] [ j ] ) ;

} // оценка квадрата нормы резольвентного оператора R s t d : : c o u t "Оцениваемквадратнормырезольвентногооператора... " ;

std : : cout. f l u s h ( ) ;

MATH: : INTERVAL sum ( 0. 0 ) ;

f o r ( s i z e _ t j 1 = 0 ;

j 1 N;

j 1++) f o r ( s i z e _ t j 2 = 0 ;

j 2 N;

j 2++) f o r ( s i z e _ t i = 0 ;

i N;

i ++) f o r ( s i z e _ t k = 0 ;

k N;

k++) { MATH: : INTERVAL d = invA_int ( i, j 1 ) invA_int ( i, j 2 ) ;

d = d K _ coeff_ int ( j 1, k ) K _ coeff_ int ( j 2, k ) ;

d = d norm2_int [ c o e f f _ r d [ i ]. d e g r e e ] ;

d = d norm2_int [ c o e f f _ r d [ k ]. d e g r e e ] ;

d = d ( tau_ int [ c o e f f _ r d [ i ]. r e g i o n +1] tau_ int [ c o e f f _ r d [ i ]. r e g i o n ] ) ;

d = d ( tau_ int [ c o e f f _ r d [ k ]. r e g i o n +1] tau_ int [ c o e f f _ r d [ k ]. r e g i o n ] ) ;

sum = sum + d ;

} mpq_class denom = ( b a s i s _ f n [0] h i ( ) b a s i s _ f n [0] l o ( ) ) ;

denom = denom ;

sum = sum / MATH: : INTERVAL( denom ) ;

s t d : : c o u t "Ok\n" ;

return sum ;

} //================================== MATH: : DOUBLE Z s o l u t i o n ( s t d : : v a l a r r a y MATH: : DOUBLE& t ) { size_t region = ;

i f ( t [ 0 ] = tau [ n ]. get_d ( ) ) r e g i o n = n1;

else { f o r ( s i z e _ t i = 1 ;

i = n ;

i ++) { i f ( t [ 0 ] tau [ i ]. get_d ( ) ) { r e g i o n = i 1;

break ;

} } } return z [ r e g i o n ] (MATH: : p h i ( t [ 0 ], MATH: : DOUBLE( tau [ r e g i o n ]. get_d ( ) ), MATH: : DOUBLE( tau [ r e g i o n + 1 ]. get_d ( ) ), b a s i s _ f n _ d b l [0] l o ( ), b a s i s _ f n _ d b l[0] h i ( ) ) ) ;

} MATH: : DOUBLE X s o l u t i o n ( s t d : : v a l a r r a y MATH: : DOUBLE& t ) { size_t region = 0;

i f ( t [ 0 ] = tau [ n ]. get_d ( ) ) r e g i o n = n1;

else { f o r ( s i z e _ t i = 1 ;

i = n ;

i ++) { i f ( t [ 0 ] tau [ i ]. get_d ( ) ) { r e g i o n = i 1;

break ;

};

} } return x [ r e g i o n ] (MATH: : p h i ( t [ 0 ], MATH: : DOUBLE( tau [ r e g i o n ]. get_d ( ) ), MATH: : DOUBLE( tau [ r e g i o n + 1 ]. get_d ( ) ), b a s i s _ f n _ d b l [0] l o ( ), b a s i s _ f n _ d b l[0] h i ( ) ) ) ;

} //================================== MATH: : DOUBLE K_approx ( s t d : : v a l a r r a y MATH: : DOUBLE & t ) { // Определяем подобласть s i z e _ t regT =0, r e g S =0;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { i f ( t [ 0 ] = tau [ i + 1 ]. get_d ( ) ) { regT = i ;

break ;

} } f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { i f ( t [ 1 ] = tau [ i + 1 ]. get_d ( ) ) { r e g S = i ;

break ;

} } // [i+1, i ] [, ] MATH: : DOUBLE T = MATH: : p h i ( t [ 0 ],MATH: : DOUBLE( tau [ regT ]. get_d ( ) ), MATH: : DOUBLE( tau [ regT + 1 ]. get_d ( ) ), b a s i s _ f n _ d b l [0] l o ( ), b a s i s _ f n _ d b l[0] h i ( ) ) ;

MATH: : DOUBLE S = MATH: : p h i ( t [ 1 ],MATH: : DOUBLE( tau [ r e g S ]. get_d ( ) ), MATH: : DOUBLE( tau [ r e g S + 1 ]. get_d ( ) ), b a s i s _ f n _ d b l [0] l o ( ), b a s i s _ f n _ d b l[0] h i ( ) ) ;

MATH: : DOUBLE K_sum = 0 ;

f o r ( s i z e _ t i = 0 ;

i c o e f f _ r d. s i z e ( ) ;

i ++) { i f ( c o e f f _ r d [ i ]. r e g i o n != regT ) continue ;

MATH: : DOUBLE bfT = ( b a s i s _ f n _ d b l [ c o e f f _ r d [ i ]. d e g r e e ] ) ( T ) ;

f o r ( s i z e _ t j = 0 ;

j c o e f f _ r d. s i z e ( ) ;

j ++) { i f ( c o e f f _ r d [ j ]. r e g i o n != r e g S ) continue ;

MATH: : DOUBLE c o e f f = K_coeff [ i ] [ j ]. get_d ( ) ;

K_sum += bfT ( b a s i s _ f n _ d b l [ c o e f f _ r d [ j ]. d e g r e e ] ) ( S ) c o e f f ;

} } return K_sum;

} MATH: : DOUBLE f_approx ( s t d : : v a l a r r a y MATH: : DOUBLE & t ) { // Определяем подынтервал size_t reg = 0;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { i f ( t [ 0 ] = tau [ i + 1 ]. get_d ( ) ) { r e g = i ;

break ;

} } // [i+1, i ] [, ] MATH: : DOUBLE T = MATH: : p h i ( t [ 0 ],MATH: : DOUBLE( tau [ r e g ]. get_d ( ) ), MATH: : DOUBLE( tau [ r e g + 1 ]. get_d ( ) ), b a s i s _ f n _ d b l [0] l o ( ), b a s i s _ f n _ d b l[0] h i ( ) ) ;

MATH: : DOUBLE f_sum = 0 ;

f o r ( s i z e _ t i = 0 ;

i c o e f f _ r d. s i z e ( ) ;

i ++) { i f ( c o e f f _ r d [ i ]. r e g i o n != r e g ) continue ;

mpq_class c o e f f = f _ c o e f f [ 0 ] [ i ] ;

f o r ( s i z e _ t l = 1 ;

l = lambda. s i z e ( ) ;

l ++) { c o e f f += f _ c o e f f [ l ] [ i ] lambda [ l 1 ] ;

} f_sum += ( b a s i s _ f n _ d b l [ c o e f f _ r d [ i ]. d e g r e e ] ) ( T) c o e f f. get_d ( ) ;

} return f_sum ;

} //================================== void ub_spectrum ( mpq_class X0) { size_t N = coeff_rd. s i z e ( ) ;

// A = E A MATH: : MATRIX mpq_class t i l d e A (N,N ) ;

f o r ( s i z e _ t i = 0 ;

i N;

i ++) f o r ( s i z e _ t j = 0 ;

j N;

j ++) { t i l d e A ( i, j ) = ( i==j ? 1 : 0 ) A( i, j ) ;


} s t d : : c o u t "ВычислениехарактеристическогомногочленадляEA... " ;

std : : cout. f l u s h ( ) ;

MATH: : POLYNOM mpq_class p = t i l d e A. c h a r p o l y ( ) ;

MATH: : POLYNOM mpq_class dp = d i f f ( p ) ;

s t d : : c o u t "Ok\n" ;

mpq_class sigma_max = X0 ;

s t d : : c o u t "Всекорнивнутрикругарадиуса " X0. get_d ( ) s t d : : e n d l ;

s t d : : c o u t "Оцениваемверхнююграницуспектра... " ;

std : : cout. f l u s h ( ) ;

mpq_class e r r ( 1 e 5);

mpq_class sigma_prev = sigma_max ;

while ( true ) { sigma_max = p ( sigma_max ) / dp ( sigma_max ) ;

mpq_class s t e p = sigma_prevsigma_max ;

i f (MATH: : ABS( s t e p ) e r r ) break ;

sigma_prev = sigma_max ;

// уменьшаем размер числителя и знаменателя, // используя округление сверху MATH: : INTERVAL sigma_ int ( sigma_max ) ;

sigma_max = mpq_class ( double ( sigma_ int. h i ( ) ) ) ;

} s t d : : c o u t "Ok\n" ;

s t d : : c o u t "max sigma = " sigma_max. get_d ( ) s t d : : e n d l ;

} };

namespace PROBLEM { //================================== extern mpq_class alpha1, a l p h a 2 ;

mpq_class normZ2, normK2 ;

// Приближённые решения std : : vector MATH: : POLYNOM mpq_class z_polynom, x_polynom ;

void s o l v e _ L a g r a n g e ( ) { s t d : : c o u t "ВычислениемножителейЛагранжа... " ;

std : : cout. f l u s h ( ) ;

// интеграл по всему отрезку s t d : : v a l a r r a y mpq_class S ( mpq_class ( 0 ), lambda. s i z e ( ) + 1 ) ;

mpq_class &o r t h o = b a s i s _ f n [ 0 ] ;

MATH: : ORTHO f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { mpq_class d = ( tau [ i +1] tau [ i ] ) / ( o r t h o. h i ( ) o r t h o. l o ( ) ) ;

f o r ( s i z e _ t L = 0 ;

L = lambda. s i z e ( ) ;

L++) { MATH: : POLYNOM mpq_class p = i n t e g r a l ( z_polynom_lambda ( i, L ) ) ;

S [ L ] += ( p ( o r t h o. h i ()) p ( o r t h o. l o ( ) ) ) d ;

} } // Вычисляем множитель Лагранжа lambda [ 0 ] = ( alpha2a l p h a 1 S [ 0 ] ) / S [ 1 ] ;

s t d : : c o u t "Ok\n" ;

s t d : : c o u t " lambda = " MATH: : INTERVAL( lambda [ 0 ] ) s t d : : e n d l ;

// Собираем z приближённое решение интегрального уравнения s t d : : c o u t "Собираемприближённоерешениеинтегральногоуравнения... " ;

std : : cout. f l u s h ( ) ;

z_polynom. r e s i z e ( n ) ;

z. r e s i z e (n ) ;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { z_polynom [ i ] = z_polynom_lambda ( i, 0 ) ;

f o r ( s i z e _ t L = 1 ;

L = lambda. s i z e ( ) ;

L++) { z_polynom [ i ] += lambda [ L1] z_polynom_lambda ( i, L ) ;

} z [ i ]. r e s i z e ( z_polynom [ i ]. d e g r e e ( ) ) ;

f o r ( s i z e _ t j = 0 ;

j = z [ i ]. d e g r e e ( ) ;

j ++) { z [ i ] [ j ] = z_polynom [ i ] [ j ]. get_d ( ) ;

} } s t d : : c o u t "Ok\n" ;

// Оцениваем норму z mpq_class normZ2 = 0 ;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { MATH: : POLYNOM mpq_class p = MATH: : i n t e g r a l ( z_polynom [ i ] z_polynom [ i ] ) ;

mpq_class s = p ( o r t h o. h i ()) p ( o r t h o. l o ( ) ) ;

s = ( tau [ i +1] tau [ i ] ) ;

normZ2 += s ;

} normZ2 /= o r t h o. h i () o r t h o. l o ( ) ;

s t d : : c o u t " | | t i l d e z | | = " MATH: : SQRT(MATH: : INTERVAL( normZ2 ) ) s t d : : e n d l ;

// Собираем приближённое решение вариационной задачи t // x(t) = 1 + a z(s) ds s t d : : c o u t "Собираемприближённоерешениевариационнойзадачи... " ;

std : : cout. f l u s h ( ) ;

x_polynom. r e s i z e ( n ) ;

x. r e s i z e (n ) ;

MATH: : POLYNOM mpq_class x0 ;

x0 [ 0 ] = a l p h a 1 ;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { mpq_class d = ( tau [ i +1] tau [ i ] ) / ( o r t h o. h i () o r t h o. l o ( ) ) ;

x_polynom [ i ] = i n t e g r a l ( z_polynom [ i ] ) d ;

x0 [ 0 ] = x_polynom [ i ] ( o r t h o. l o ( ) ) ;

x_polynom [ i ] += x0 ;

x [ i ]. r e s i z e ( x_polynom [ i ]. d e g r e e ( ) ) ;

f o r ( s i z e _ t j = 0 ;

j = x [ i ]. d e g r e e ( ) ;

j ++) { x [ i ] [ j ] = x_polynom [ i ] [ j ]. get_d ( ) ;

} x0 [ 0 ] = x_polynom [ i ] ( o r t h o. h i ( ) ) ;

} s t d : : c o u t "Ok\n" ;

// Оцениваем норму x mpq_class Xsum = 0 ;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { MATH: : POLYNOM mpq_class p = MATH: : i n t e g r a l ( x_polynom [ i ] x_polynom [ i ] ) ;

mpq_class s = p ( o r t h o. h i ()) p ( o r t h o. l o ( ) ) ;

s = ( tau [ i +1] tau [ i ] ) ;

Xsum += s ;

} Xsum /= o r t h o. h i () o r t h o. l o ( ) ;

s t d : : c o u t " | | t i l d e x | | = " MATH: : SQRT(MATH: : INTERVAL(Xsum ) ) s t d : : e n d l ;

// Оценка норм K и f s t d : : c o u t "Оценканорм приближений... " ;

std : : cout. f l u s h ( ) ;

// интервал, на котором // базисные функции ортогональны mpq_class BA = b a s i s _ f n [0] h i () b a s i s _ f n [0] l o ( ) ;

// Оценка f mpq_class normF2 = 0 ;

f o r ( s i z e _ t i = 0 ;

i c o e f f _ r d. s i z e ( ) ;

i ++) { mpq_class f c = f _ c o e f f [ 0 ] [ i ] ;

f o r ( s i z e _ t L = 1 ;

L = lambda. s i z e ( ) ;

L++) { f c += f _ c o e f f [ L ] [ i ] lambda [ L 1 ] ;

} fc = fc fc ;

f c = tau [ c o e f f _ r d [ i ]. r e g i o n +1] tau [ c o e f f _ r d [ i ]. r e g i o n ] ;

f c = b a s i s _ f n [ c o e f f _ r d [ i ]. d e g r e e ]norm2 ( ) ;

normF2 += f c ;

} normF2 /= BA;

// Оценка K normK2 = 0 ;

for ( size_t i = 0 ;

i coeff_rd. s i z e ( ) ;

i ++) { mpq_class SK = 0 ;

f o r ( s i z e _ t j = 0 ;

j c o e f f _ r d. s i z e ( ) ;

j ++) { mpq_class Kc = K_coeff [ i ] [ j ] K_coeff [ i ] [ j ] ;

Kc = tau [ c o e f f _ r d [ j ]. r e g i o n +1] tau [ c o e f f _ r d [ j ]. r e g i o n ] ;

Kc = b a s i s _ f n [ c o e f f _ r d [ j ]. d e g r e e ]norm2 ( ) ;

SK += Kc ;

} SK = tau [ c o e f f _ r d [ i ]. r e g i o n +1] tau [ c o e f f _ r d [ i ]. r e g i o n ] ;

SK = b a s i s _ f n [ c o e f f _ r d [ i ]. d e g r e e ]norm2 ( ) ;


normK2 += SK ;

} normK2 /= (BABA ) ;

s t d : : c o u t "Ok\n" ;

s t d : : c o u t " | | t i l d e F | | = " MATH: : SQRT(MATH: : INTERVAL( normF2 ) ) s t d : : e n d l ;

s t d : : c o u t " | | t i l d e K | | = " MATH: : SQRT(MATH: : INTERVAL( normK2 ) ) s t d : : e n d l ;

} };

main.cpp #include i o s t r e a m #include f s t r e a m #include time. h #include " problem. hpp" i n t main ( ) { using namespace s t d ;

using namespace PROBLEM;

cout. p r e c i s i o n ( 2 0 ) ;

clock_t s t a r t = c l o c k ( ) ;

try { init ();

//============================================= // аппроксимируем ядро и правую часть create_coeffs ( ) ;

// Коррекция Kij f o r ( s i z e _ t i = 0 ;

i c o e f f _ r d. s i z e ( ) ;

i ++) f o r ( s i z e _ t j = 0 ;

j c o e f f _ r d. s i z e ( ) ;

j ++) { mpq_class s = ( K_coeff [ i ] [ j ]+ K_coeff [ j ] [ i ] ) / 2 ;

K_coeff [ i ] [ j ] = K_coeff [ j ] [ i ] = s ;

} // строим СЛАУ create_matrix_A ( ) ;

create_ v ector_ B ( ) ;

// решаем СЛАУ и находим множители Лагранжа solve_IEF ( ) ;

solve_Lagrange ( ) ;

//============================================= // генерируем данные для графиков c o u t "Генерируемданныедляграфиков... " ;

cout. f l u s h ( ) ;

const s i z e _ t graphN = 1 0 ;

valarray MATH: : DOUBLE t ;

// f, f, f f, z, x o f s t r e a m ofs_ t, o f s _ f, ofs_f_approx, ofs_ df, ofs_z, ofs_x ;

o f s _ t. open ( ( name+"_t" ). c _ s t r ( ) ) ;

o f s _ f. open ( ( name+"_f" ). c _ s t r ( ) ) ;

ofs_ f_ approx. open ( ( name+" _f_approx" ). c _ s t r ( ) ) ;

o f s _ d f. open ( ( name+"_df" ). c _ s t r ( ) ) ;

ofs_ z. open ( ( name+"_z" ). c _ s t r ( ) ) ;

ofs_x. open ( ( name+"_x" ). c _ s t r ( ) ) ;

t. resize (1);

o f s _ t ( n ( graphN +1)) e n d l ;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { MATH: : DOUBLE a = tau [ i ]. get_d ( ), b = tau [ i + 1 ]. get_d ( ) ;

MATH: : DOUBLE h = ( ba ) / graphN ;

f o r ( s i z e _ t j = 0 ;

j = graphN ;

j ++) { t [ 0 ] = a + h j ;

F_LAMBDA lambda_num ( 0 ) ;

MATH: : DOUBLE f _ v a l = f ( t,&lambda_num ) ;

f o r ( s i z e _ t l = 1 ;

l = lambda. s i z e ( ) ;

l ++) { lambda_num. number = l ;

f _ v a l += f ( t,&lambda_num ) lambda [ l 1 ]. get_d ( ) ;

} MATH: : DOUBLE f_approx_val = f_approx ( t ) ;

o f s _ t t [ 0 ] e n d l ;

o f s _ f f _ v a l e n d l ;

ofs_ f_ approx f_approx_val e n d l ;

o f s _ d f ( f_ v alf_approx_val ) e n d l ;

ofs_ z Z s o l u t i o n ( t ) e n d l ;

ofs_x X s o l u t i o n ( t ) e n d l ;

} } ofs_x. c l o s e ( ) ;

ofs_ z. c l o s e ( ) ;

ofs_df. c l o s e ( ) ;

ofs_ f_ approx. c l o s e ( ) ;

ofs_f. clos e ( ) ;

ofs_t. c l o s e ( ) ;

// K, K, K K o f s t r e a m ofs_K, ofs_K_approx, ofs_dK ;

ofs_K. open ( ( name+"_K" ). c _ s t r ( ) ) ;

ofs_K_approx. open ( ( name+"_K_approx" ). c _ s t r ( ) ) ;

ofs_dK. open ( ( name+"_dK" ). c _ s t r ( ) ) ;

t. resize (2);

f o r ( s i z e _ t i 1 = 0 ;

i 1 n ;

i 1 ++) { MATH: : DOUBLE a t = tau [ i 1 ]. get_d ( ), bt = tau [ i 1 + 1 ]. get_d ( ) ;

MATH: : DOUBLE ht = ( bta t ) / graphN ;

f o r ( s i z e _ t j 1 = 0 ;

j 1 = graphN ;

j 1++) { t [ 0 ] = a t + ht j 1 ;

f o r ( s i z e _ t i 2 = 0 ;

i 2 n ;

i 2 ++) { MATH: : DOUBLE a s = tau [ i 2 ]. get_d ( ), bs = tau [ i 2 + 1 ]. get_d ( ) ;

MATH: : DOUBLE hs = ( bsa s ) / graphN ;

f o r ( s i z e _ t j 2 = 0 ;

j 2 = graphN ;

j 2++) { t [ 1 ] = a s + hs j 2 ;

MATH: : DOUBLE K_val = K( t ) ;

MATH: : DOUBLE K_approx_val = K_approx ( t ) ;

ofs_K K_val e n d l ;

ofs_K_approx K_approx_val e n d l ;

ofs_dK ( K_valK_approx_val) e n d l ;

} // j } // i } // j } // i ofs_dK. c l o s e ( ) ;

ofs_K_approx. c l o s e ( ) ;

ofs_K. c l o s e ( ) ;

c o u t "Ok\n" ;

//============================================= // оцениваем (грубо!) погрешность аппроксимации const s i z e _ t errN = 1 0 0 ;

// погрешность f f c o u t "Оцениваемпогрешностьаппроксимации f... " ;

t. resize (1);

std : : valarray MATH: : DOUBLE DeltaF ( n ) ;

MATH: : DOUBLE normDeltaF = 0 ;

f o r ( s i z e _ t i = 0 ;

i n ;

i ++) { DeltaF [ i ] = 0 ;

MATH: : DOUBLE a = tau [ i ]. get_d ( ), b = tau [ i + 1 ]. get_d ( ) ;

MATH: : DOUBLE h = ( ba ) / errN ;

f o r ( s i z e _ t j = 0 ;

j = errN ;

j ++) { t [ 0 ] = a + h j ;

F_LAMBDA lambda_num ( 0 ) ;

MATH: : DOUBLE d f = f ( t,&lambda_num ) ;

f o r ( s i z e _ t l = 1 ;

l = lambda. s i z e ( ) ;

l ++) { lambda_num. number = l ;

d f += f ( t,&lambda_num ) lambda [ l 1 ]. get_d ( ) ;

} d f = f_approx ( t ) ;

i f ( DeltaF [ i ] MATH: : ABS( d f ) ) DeltaF [ i ] = MATH: : ABS( d f ) ;

} normDeltaF += DeltaF [ i ] DeltaF [ i ] ( ba ) ;

} normDeltaF = s q r t ( normDeltaF ) ;

c o u t "Ok\n" ;

c o u t "max DeltaF =" DeltaF. max ( ) e n d l ;

c o u t " | | DeltaF | | =" normDeltaF e n d l ;

// погрешность K K c o u t "ОцениваемпогрешностьаппроксимацииK... " ;

t. resize (2);

MATH: : MATRIX MATH: : DOUBLE DeltaK ( n, n ) ;

MATH: : DOUBLE maxDeltaK = 0 ;

MATH: : DOUBLE normDeltaK = 0 ;

f o r ( s i z e _ t i 1 = 0 ;

i 1 n ;

i 1 ++) f o r ( s i z e _ t i 2 = 0 ;

i 2 n ;

i 2 ++) { DeltaK ( i 1, i 2 ) = 0 ;

MATH: : DOUBLE a t = tau [ i 1 ]. get_d ( ), bt = tau [ i 1 + 1 ]. get_d ( ) ;

MATH: : DOUBLE ht = ( bta t ) / errN ;

MATH: : DOUBLE a s = tau [ i 2 ]. get_d ( ), bs = tau [ i 2 + 1 ]. get_d ( ) ;

MATH: : DOUBLE hs = ( bta t ) / errN ;

f o r ( s i z e _ t j 1 = 0 ;

j 1 = errN ;

j 1++) f o r ( s i z e _ t j 2 = 0 ;

j 2 = errN ;

j 2++) { t [ 0 ] = a t + ht j 1 ;

t [ 1 ] = a s + hs j 2 ;

MATH: : DOUBLE dK = MATH: : ABS(K( t )K_approx ( t ) ) ;

i f ( DeltaK ( i 1, i 2 ) dK) DeltaK ( i 1, i 2 ) = dK ;

} i f ( maxDeltaK DeltaK ( i 1, i 2 ) ) maxDeltaK = DeltaK ( i 1, i 2 ) ;

normDeltaK += DeltaK ( i 1, i 2 ) DeltaK ( i 1, i 2 ) ( bta t ) ( bsa s ) ;

} normDeltaK = s q r t ( normDeltaK ) ;

c o u t "Ok\n" ;

c o u t "max DeltaK =" maxDeltaK e n d l ;

c o u t " | | DeltaK | | =" normDeltaK e n d l ;

//============================================= // оцениваем норму резольвентного оператора...

MATH: : INTERVAL normR = SQRT( estimate_N2RC ( ) ) ;

c o u t " | | R | | = " normR e n d l ;

//...и проверяем выполнение теоремы об обратном операторе MATH: : INTERVAL one ( 1. 0 ) ;

c o u t " 1 / ( 1 + | |R | | ) = " ( one / ( one+normR ) ) e n d l ;

s t d : : c o u t " | | z t i l d e z | | =" (1+normR. h i ())/(1 (1+ normR. h i ( ) ) normDeltaK ) ( normDeltaK s q r t ( normZ2. get_d ())+ normDeltaF ) s t d : : e n d l ;

//============================================= // Оцениваем верхнюю границу спектра K MATH: : DOUBLE X0 = MATH: : SQRT(MATH: : INTERVAL( normK2 ) ). h i ( ) + normDeltaK 2 ;

ub_spectrum ( mpq_class ( double (X0 ) ) ) ;

// */ destroy ( ) ;

c o u t "Успешноезавершение\n" ;

} catch ( const char msg ) { c o u t "Ошибка : " msg e n d l ;

} cout. p r e c i s i o n ( 5 ) ;

c o u t "Затрачено " f l o a t ( c l o c k () s t a r t ) /CLOCKS_PER_SEC "сек. \ n" ;

return 0 ;

} Приложение B GNU Octave. Вспомогательные процедуры B.1 Решение скалярного уравнения Метод деления отрезка пополам (медленный, но надёжный):

function r e s = NLEsolver ( l, r, f, eps = 1 e 5) f l = feval ( f, l ) ;

f r = feval ( f, r ) ;

i f ( f l == 0 ) r e s = l ;

return ;

e ndif i f ( f r == 0 ) r e s = r ;

return ;

e ndif ;

i f ( f l f r 0) error ( " NLEsolver : f ( l ) f ( r ) 0 " ) ;

e ndif while ( rl eps ) m = ( l+r ) / 2 ;

fm = f e v a l ( f,m) ;

i f ( fm == 0 ) r e s = m;

return ;

e ndif i f ( f l fm 0 ) l = m;

f l = fm ;

else r = m;

f r = fm ;

e ndif endwhile r e s = ( l+r ) / 2 ;

endfunction B.2 Решение задачи Коши Решение задачи Коши x (t) = f t, x(t), x(a) = с помощью одношаговых методов в общем случае может быть представлено в следующем виде:

x0 =, xi+1 = xi + h(ti, xi ), Для метода Рунге–Кутта 4-го порядка (t, x) = (K1 + 2K2 + 2K3 + K4 ), где K1 =f (t, x), K2 =f (1 + h/2, x + hK1 /2), K3 =f (1 + h/2, x + hK2 /2), K4 =f (1 + h, x + hK3 ) function r e s = RK4 ( t, x, f, h ) k1 = f e v a l ( f, t, x ) ;

k2 = f e v a l ( f, t+h / 2, x+h /2 k1 ) ;

k3 = f e v a l ( f, t+h / 2, x+h /2 k2 ) ;

k4 = f e v a l ( f, t+h, x+h k3 ) ;

r e s = 1 / 6 ( k1 +2( k2+k3)+k4 ) ;

endfunction Применение одношагового метода:

function r e s = ODEsolver ( t0, tn, x0, f, n, method = "RK4" ) h = ( tnt 0 ) / n ;

x = zeros ( n, length ( x0 ) ) ;

x ( 1, : ) = x0 ;

f o r i = 1 : ( n1) x ( i + 1, : ) = x ( i, : ) + h f e v a l ( method, t 0+h i, x ( i, : ), f, h ) ;

endfor res = x ;

endfunction B.3 Решение краевой двухточечной задачи d2 x(t) 1 + cos(t) = sin(2t) exp(t/5)x(t) +, dt2 x(0) = 1, x(10) = Правая часть системы обыкновенных дифференциальных уравнений пер вого порядка:

function r e s = f 2 ( t, x ) r e s = [ x ( 2 ), sin ( 2 t ) exp( t / 5 ) x (1)+(1+ cos ( t ) ) / 2 ] ;

endfunction Функция, возвращающая значение x(10) в зависимости от x (0) = 2 :

function r e s = ODE ( a l p h a ) x = ODEsolver ( 0, 1 0, [ 1, a l p h a ], " f 2 ", 1 0 0 0 ) ;

re s = x( size (x )(1),1) 1;

endfunction Решение нелинейного уравнения x 10 | 2 = и построение графика x(t) для найденного значения 2 :

a l p h a = NLEsolver ( 5,1, "ODE" ) x = ODEsolver ( 0, 1 0, [ 1, a l p h a ], " f 2 ", 1 0 0 0 ) ;

alpha = -2. plot ( ( 1 : length ( x ( :, 1 ) ) ) 0. 0 1, x ( :, 1 ) ’, "b", " l i n e w i d t h ", 2 ) grid on print d e p s c 2 oct 03. eps B.4 Графики B.4.1 2D problem = "ex_1_15_3" ;

g r = [ "_f" ;

"_df" ;

" _f_approx" ;

"_z" ;

"_x" ] ;

fid_ x = fopen ( s t r c a t ( ".. / prog / data / ", problem, "_t" ), " r " ) ;

n = f s c a n f ( fid_x, "%d", 1 ) ;

x = zeros ( 1, n ) ;

f o r i =1: n x ( i ) = f s c a n f ( fid_x, "%f ", 1 ) ;

endfor f c l o s e ( fid_ x ) ;

for g = 1 : ( size ( gr ) ( 1 ) ) fid_ y = fopen ( s t r c a t ( ".. / prog / data / ", problem, deblank ( g r ( g, : ) ) ), " r " ) ;

y = zeros ( 1, n ) ;

f o r i =1: n y ( i ) = f s c a n f ( fid_y, "%f ", 1 ) ;

endfor f c l o s e ( fid_ y ) ;

plot ( x, y, "b", " l i n e w i d t h ", 2 ) grid on f i l e n a m e = s t r c a t ( problem, deblank ( g r ( g, : ) ), ". e p s " ) ;

print ( f i l e n a m e, ’d e p s c 2 ’ ) endfor B.4.2 3D problem = "ex_1_15_3" ;

g r = [ "_K" ;

"_dK" ;

"_K_approx" ] ;

fid_ x = fopen ( s t r c a t ( ".. / prog / data / ", problem, "_t" ), " r " ) ;

n = f s c a n f ( fid_x, "%d", 1 ) ;

x = zeros ( 1, n ) ;

f o r i =1: n x ( i ) = f s c a n f ( fid_x, "%f ", 1 ) ;

endfor f c l o s e ( fid_ x ) ;

X= zeros ( n, n ) ;

i =1: n X ( :, i )=x ;

endfor for Y= zeros ( n, n ) ;

i =1: n Y( i, : ) = x ;

endfor for Z = zeros ( n, n ) ;

f o r g = 1:( s i z e ( g r ) ( 1 ) ) f i d _ z = fopen ( s t r c a t ( ".. / prog / data / ", problem, deblank ( g r ( g, : ) ) ), " r " ) ;

f o r i =1: n f o r j =1:n Z ( i, j ) = f s c a n f ( fid_ z, "%f ", 1 ) ;

endfor endfor fclo se ( fid_z ) ;

surf (X, Y, Z ) view ( 4 0, 3 0 ) f i l e n a m e = s t r c a t ( problem, deblank ( g r ( g, : ) ), ". e p s " ) ;

print ( f i l e n a m e, ’d e p s c 2 ’ ) c o n t o u r f ( x, x, Z) f i l e n a m e = s t r c a t ( problem, deblank ( g r ( g, : ) ), "_im. e p s " ) ;

print ( f i l e n a m e, ’d e p s c 2 ’ ) endfor

Pages:     | 1 | 2 ||
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.