Курсовая работа по дисциплине Дифференциальные уравнения

 Федеральное государственное бюджетноеобразовательное учреждениевысшего профессионального образования«Московский государственный технический университетим. Н. Э. Баумана».Курсовая работапо предмету «Дифференциальные уравнения»Выполнил:студент группы АК3-31аэрокосмического факультетаСкрябин ДмитрийПреподаватель: Апельцин В.Ф.IПостановка задачи.1. Математический аспектДанная курсовая работа посвящена построению приближенного решения краевой задачидля дифференциального уравнения вида:d2 u+n( x )u ( x )=0dx 2По условию, n(x)представляет собой кусочно-непрерывную функцию, равную константеk2 на интервалах (-∞; 0) и (а; +∞) и непрерывной известной функции на [0; a].

На границах[0; a] решение и его первая производная должны удовлетворять условиям:(гладкое сшивание решения)u(-0) = u(+0)u(a-0) = u(a+0)u’(-0) = u’(+0)u’(a-0) = u’(a+0)2. Физический аспектЗадача соответствует возбуждению плоского слоя неоднородной намагниченной средыэлектромагнитным полем, вектор электрического поля распространяется вдоль оси х.

kимеет смысл волнового числа в однородной части пространства, n( x ) - показательпреломления неоднородной среды.Из физических соображений, решение справа от слоя [0; a] должно иметь вид u( x ) =−ikxBe,ikx−ikxа слева u( x ) = e +Ae. Задача состоит в приближенном построении решениявнутри неоднородного слоя и нахождении комплексных констант A и B (называемых,соответственно, коэффициентами отражения и прохождения электромагнитного поля).n( x )Действительная область значений функциисоответствует отсутствиюпоглощения энергии в среде, так что критерием правильности полученных приближенныхзначений служит энергетическое тождество:|A|2 +|B|2 =11.

Постановка задачи.Распространяющаяся в однороднойдифференциальным уравнением видаd2 y 2+k y= 0dx 2изотропнойсредеВнутрь среды помещают неоднородный слой ширинойуравнение распространения волны удовлетворяет уравнению2d y+n( x ) y= 0dx 2Требуется найти функциюпросачивания волны.y( x)волнаописывается(1)a . Внутри этого слоя(2)внутри слоя и коэффициенты отражения и2. Теоретическая часть.OxВведём систему координат таким образом, что осьсовпадает сOyнаправлением распространения волны, осьей перпендикулярна, а прямая x=0совпадает с левой границей слоя. Тогда для правой границы x=a .

В этом случае можноk 2 ,x ∉ ( 0, a )n ( x ) ,0 <x<aN ( x ) =¿ {}N ( x ) , для которойввести функциюбудем рассматривать уравнениеd2 u+N ( x ) u=0dx 2.Тогдаikx(3)К слою подходит волна, описываемая функцией Y ( x )=e, являющейсярешением уравнения (1). На границе слоя часть волны отразится. С учётом изотропностиx<0слоя можно записать, что прирешением уравнения (3) будет функцияikx−ikxY ( x )=e +Ae, где A – числовой коэффициент, называемый коэффициентом−ikxотражения. Аналогично, при x>a решение примет вид Y ( x )=Be, где B –числовой коэффициент, называемый коэффициентом просачивания.Допустим, нам удалось получить решение уравнения (3) на [0 ;a ] .

По теореме оструктуре решения обыкновенного линейного дифференциального уравнении, этоY ( x )=c 1 y 1 ( x )+c 2 y 2 ( x ) . На решение наложим условиерешение имеет виднепрерывности и гладкости в точках x=0 иY (0 )|−0 =Y ( 0 )|+0Y' (0 )|−0 =Y' ( 0 )|+0Y (a )|−0 =Y ( a )|+0Y' (a )|−0 =Y' ( a)|+0x=a . Они имеют вид{}{}{}(4)Систему (4) можно привести, с учётом известного нам поведения функцииU ( x ) , к видуY ( 0 )=1 +AY'( 0 )=ik ( 1− A )Y (a )=Be ikaY' ( a) =ikBeika{}{}{}(5)Введём функциюz ( x )=dydx .

Тогда уравнение (2) сводится к системеdz+n ( x ) y= 0dxdy=zdx{}(6)Систему (6) представим в виде векторного уравненияddxy (x)z(x)ri gh(0( −n()=y (z(x)10)x)x)ri gh()(7)Символически запишем (7) в видеd ⃗⃗W =D ( x )⋅Wdx(8)Введём фундаментальную матрицу Ω( x ) , удовлетворяющую условиямd Ω( x )=D( x )⋅Ω ( x )dxΩ ( 0) =E{}(9)Тогда очевидно, что⃗ ( x )=Ω ( x )⋅W⃗ (0 )W(10)Из (5) видно, что1 +Aik ( 1 − A )righ( )( 0 )=⃗W(11)Раскроем выражение (10) с учётом (11):y(xz (x))ri ghik1 +A( 1 − A)ri gh()(12)Подставимx=aв (12):yz(a(a))rig hik1 +A( 1 −A)rig hω 11 (ω21(11+A+A))+ ω12+ ω 22ikik(1 −(1−AA))rig h()(13)В то же время, в силу (5) имеем:y(z(a )a )rig hB eik aikB e i karig h()(14)Тогда из (13) и (14) получим систему из двух линейных уравнений скомплексными коэффициентами относительно A и B :( ω 11 ( a )−ikω12 ( a )) A −eika B=−ω11 ( a )−ikω12 ( a )( ω21 ( a )−ikω 22 ( a) ) A −ike ika B=− ω21 ( a)−ikω22 ( a ){}(15)Таким образом, для определения коэффициентов отражения и просачивания намдостаточно найти значение фундаментальной матрицы при x=a , т.е.

в крайней точкеинтервала интегрирования системы (2). Построение явного решения Y ( x ) потребуетпостроения явного вида фундаментальной матрицы. Это не представляется возможным,поэтому мы воспользуемся численным методом для вычисления приближённых значенийΩ( x ) во внутренних точках отрезка x ∈[ 0 ;a ] .3. Построение численного решения.Разобьём отрезок x ∈[ 0 ;a ] на M равных отрезков вида [ x j ;x j+1 ] , гдеax j= jM, j= 0, M .

На каждом из отрезков аппроксимируем функцию n( x )a1a~x j = ( j+ )h=~~n(x)=n(x)M2 . Примем шагM , тогдаjпостоянной, гдеh~n ( x )=h⋅j+2 . Тогда на [ x j ;x j+1 ] матрица A ( x ) в уравнении (8) будет неhA j ( x ) =A( x j + )2 . Тогда на [ 0 ;h]переменной, а постоянной. Будем считатьD0 hD1 h D0 h[0 ; 2 h ]ee ⋅eрешением будет, наи т.д. Таким образом,фундаментальная матрица может быть найдена по формуле.MΩ(a )=∏ e A hj(16)После этого из системы (15) определяются коэффициенты A и B . Дляанализа корректности полученных результатов воспользуемся законом сохранения22энергии в форме |A| +|B| =1 . При численном решении строгое равенство перейдёт вj= 1приближённое.

Для построения функции Y ( x ) нужно найти её значения в узловыхточках, а затем выполнить интерполяцию. Для нахождения значения в точке x= ¿ ¿воспользуемся указанным выше алгоритмом, но не на отрезке [0 ;a ] , а на отрезке[ 0, ¿ ¿ .4. Исходные данные.При решении конкретной задачи были использованы следующие исходныеданные:k 5, a 2 25 , x 0; 25  3 x , 0<x<1;n( x ) 2 4( x  1.5)  27 , 1 x<2; 25 , x 2.5. Реализация метода.Для численного решения задачи и наглядного представления результатов на языкес++ была написана программа, в которой реализованы класс матрица,комплексные числаи операции над ними. Сначала вычисляются значения Ω( a ) ,A,Bи точек дляпостроения графика, а затем строятся вl графики мнимой и действительной частейY ( x ) .

Полный текст программы см. в приложении.Функция нахождения матричной экспоненты:m2x2 m2x2::matr_exp(double x1, double x2){m2x2 result,E,temp(*this);double delta=x2-x1;double n2=1.0/2,n3=1.0/6,n4=1.0/24;result=E+temp*delta+temp*temp*delta*delta*n2+temp*temp*temp*delta*delta*delta*n3+temp*temp*temp*temp*delta*delta*delta*delta*n4;return result;}Функция нахождения Ω( x ) :m2x2 m2x2::fundametnal_matr(double a, double b, int N){m2x2 result, temp(*this);double h=abs(a-b)/(N);double xi=h/2;for(int i=0;i<N;i++){temp.set_A(xi);result=result*temp.matr_exp(xi-h/2,xi+h/2);xi+=h;}return result;}Функция вычисления значения искомой функции в серии точек, достаточной дляпостроения графика и записи этих значений в файл:Procedure MakeGrFile(a:real;Am:comp;N:integer); //Создание файла созначениями для графикаvarx,h:real;f,f_:textfile;j:integer;oyim,oyre:real;omx:matr;g,ed,om1,om2,ed1:comp;begined.re:=1;ed.im:=0;om1.im:=0;om2.im:=0;ed1.im:=k;ed1.re:=0;assignfile(f,'dotsre.txt');assignfile(f_,'dotsim.txt');rewrite(f);rewrite(f_);h:=a/N;x:=0;for j:=0 to N-1 dobeginomx:=FindOmega(x,1000);om1.re:=omx[1,1];om2.re:=omx[1,2];g:=plusc(CMult(om1,plusc(ed,Am)),CMult(om2,CMult(ed1,minusc(ed,Am))));oyim:=g.im;oyre:=g.re;writeln(f_,oyim:5:5);writeln(f,oyre:5:5);x:=x+h;end;closefile(f);closefile(f_);end;6.

Результаты.При решении задачи с учётом условий, указанных в пункте 4, были полученыследующие результаты:При числе разбиений 600:1)A =- 0.019521 + 0.0019084iB =0.930281 + 0.366323iA2 + B 2 =12)Графическое представление решения.Y1.21.080.960.840.720.60.480.360.240.12X0-0.12-0.24-0.36-0.48-0.6-0.72-0.84-0.96-1.08-1.20.1 0.20.3 0.4 0.5 0.6 0.7 0.8 0.911.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.92При числе разбиений 800:1)A =- 0.0195206 + 0.00190847iB =0.930281 + 0.366323iA2 + B 2 =12)Графическое представление решения.Y1.21.080.960.840.720.60.480.360.240.12X0-0.12-0.24-0.36-0.48-0.6-0.72-0.84-0.96-1.08-1.20.1 0.2 0.3 0.4 0.50.6 0.7 0.8 0.911.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.92При числе разбиений 1000:A =- 0.0195205 + 0.0019085i 1)B =0.930281 + 0.366323iA2 + B 2 =12)Графическое представление решения.Y1.21.080.960.840.720.60.480.360.240.12X0-0.12-0.24-0.36-0.48-0.6-0.72-0.84-0.96-1.08-1.20.1 0.20.3 0.4 0.5 0.6 0.7 0.8 0.911.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.92II8.

Метод последовательных приближений.Известно, что решение уравненияd2 udx 2 +k2u=f(x) (19)на всей числовой прямой -   x  приходящих из бесконечности, имеет видудовлетворяющее условию отсутствия волн,∞∫ G( x, x 1 )f ( x 1 )dx 1u(x)= −∞где G(x, x1) - функция Грина: решение той же задачи с правой частьюf(x) =  ( x – x 1 ).Для уравнения (19) с постоянным коэффициентом k2 функция Грина выписывается вявном видеi k | x − x 1|e2 ikG(x,x1)=(20)Следовательно, решение уравнения (19) выписывается в явном виде:∞ i k |x − x1|u(x)=∫−∞ e 2ikf (x1 )dx 1. (21)Перепишем уравнение (1) в эквивалентной формеd 2u k 2u (k 2  n( x))u ( x).2dx(22)Тогда, согласно (21), решение уравнения (22) с правой частью (k2 - n(x))u(x)выписать в виде∞u(x)=∫ei k |x − x1|2ik−∞(k 2 - n (x 1 ))u( x 1 )dx 1.

(23)Так как (k2 - n(x))  0 при x  0, и при xзаменяется интегралом в конечных пределахau(x)=∫ei k |x − x1|2 ikможно a , то несобственный интеграл в (23)(k 2 - n ( x 1 ))u( x 1 ) dx 1.(24)Раскрывая модуль в показателе экспоненты для x  0 и для x  a, получим, что0ae−i k x∫ ei k x ( k 2 - n(x 1 ))u( x 1)dx 12ik0u(x)=для x<0,(25)иaei k x∫ e− i k x (k 2 - n ( x1 )) u( x 1 )dx 12ik0u(x)=для х>a(26)Следовательно, решение (24) содержит лишь волны, уходящие в - −iωtзависимость решения от времени выбрана в виде e).11и в (еслиikxeНо исходная постановка задачи содержит также единственную волну,приходящую из -  .ikxЕсли взять сумму полей (24) и e, то получим представление для полного поля вовсей области -   x  a∫u(x)=ei k |x − x1|2 ik0(k 2 - n ( x 1 ))u( x 1 ) dx 1eikx.(27)В сокращенных обозначениях это уравнение записываеся в видеu = Au + f.a∫e(28)i k |x − x1|2 ik(k 2 - n ( x 1 ))u( x 1 ) dx 1Здесь А – интегральный оператор 0, действующий нафункцию u(x).Нетрудно убедиться, что функция u(x), представленная в виде (27), удовлетворяет всемусловиям исходной задачи.

Действительно, правая часть (27) удовлетворяет уравнениюd2 udx 2 + k2u =0вне слоя [0, a] , и уравнению2d udx 2 + n(x)u = 0внутри этого слоя, благодаря свойствам функции Грина (20) . Краевые условия (2) такжевыполняются, что проверяется непосредственно.