Построение приближенного решения краевой задачи для дифференциального уравнения (Огромное количество решённых курсовых), страница 2

. (17)

При численном решении строгое равенство перейдёт в приближённое. Для построения функции нужно найти её значения в узловых точках, а затем выполнить интерполяцию. Для нахождения значения в точке воспользуемся указанным выше алгоритмом, но не на отрезке , а на отрезке .

7. Реализация метода

Для численного решения задачи и наглядного представления результатов на языке Delphi в программе Embarcadero Delphi XE3 Version 17.0.4625.53395 (VLC Forms Application Delphi) была написана программа, в которой реализованы матричный и комплексный типы данных и операции над ними.

8. Визуальная оболочка программы

9. Код программы c учетом исходных данных

unit Course_work_by_Yagubov_R_B_V15_AK3;

interface

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, StdCtrls, ExtCtrls,Math, Matrix, cmplx;

type

TForm1 = class(TForm)

Label5: TLabel;

Panel2 : TPanel;

PaintBox1 : TPaintBox;

PaintBox2 : TPaintBox;

Edit_Re_A : TEdit;

Edit_Re_B : TEdit;

Edit_Im_A : TEdit;

Edit_Im_B : TEdit;

Edit1: TEdit;

Edit2: TEdit;

Edit3: TEdit;

Edit4: TEdit;

Edit5: TEdit;

Edit_Check: TEdit;

procedure Edit1Change(Sender: TObject);

end;

{ объявление констант}

const

e=1e-5; //точность

interval=2; //интервал [0;interval]

k=5; //волновое число

{ объявление переменных}

var

Form1:TForm1; //форма

dx, //интервал разбиения

Min_y,Max_y,Min_y1,Max_y1,Min_y2,Max_y2:real; //мин и макc функции

N,i,j:integer; //разбиения, параметры

Coefficient_A,Coefficient_B, //коэффициент A,B

c1,c2,c3,c4,c5,c6,c7,c8:TComplex; //комплексные числа

A,matrix_eps,omega: array of MatrixPtr; //матричные массивы

u: array of TComplex; //комплексный массив

implementation

{$R *.dfm}

{ задание функции на отрезке [0,a] }

function f(x:real):real;

begin

if ( x <= 0 ) then f := 25;

if ( x > 0 ) and ( x < 2 ) then f:=3*x*sin(2*pi*x)+25;

if ( x >= 2 ) then f := 25;

end;

{результат функции u(x)}

function res_u(x:real; A:TComplex):TComplex;

begin

c1.x:=0;

c1.y:=k*x;

c1:=cExp(c1);

c2.x:=0;

c2.y:=-k*x;

c2:=cExp(c2);

c2:=cMul(A,c2);

res_u:=cSum(c1,c2);

end;

{ вычисление факториала }

function factorial(const x:integer):integer;

var

f,q: integer;

begin

f:=1;

for q:=1 to x do f:=f*q;

factorial:= f;

end;

{ вычисление разложений }

function Makloren(A:MatrixPtr; currx,x:real):MatrixPtr;

var

i,j:integer;

b,c:MatrixPtr;

begin

b:=CreateSquareMatrix(2);

SetMatrixElement(b,1,1,1);

SetMatrixElement(b,1,2,0);

SetMatrixElement(b,2,1,0);

SetMatrixElement(b,2,2,1);

for i:=1 to 4 do

begin

c:=cloneMatrix(A);

for j:=1 to i-1 do c:=MultipleMatrixOnMatrix(c,a);

c:=MultipleMatrixOnNumber(c,IntPower((currx-x),i)/factorial(i));

b:=AddMatrixOnMatrix(b,c);

end;

Makloren:=b;

end;

{ выполнение программы}

procedure TForm1.Edit1Change(Sender: TObject);

var

i,j:integer;

dx1,dy1,dx2,dy2:real;

begin

N:=StrToInt(Edit1.Text);

dx:=interval/N;

setlength(A,N);

{ построение постоянной матрицы Ai }

for i:=0 to N-1 do

begin

A[i]:=CreateSquareMatrix(2);

SetMatrixElement(A[i],1,1,0);

SetMatrixElement(A[i],1,2,1);

SetMatrixElement(A[i],2,1,-1*f((i+1)*dx-dx/2));

SetMatrixElement(A[i],2,2,0);

end;

setlength(matrix_eps,N);

{ создание массива матричных экспонент }

for i:=0 to N-1 do matrix_eps[i]:=Makloren(a[i],(i+1)*dx,i*dx);

setlength(omega,N);

{ построение массива для приближенной фундаментальной матрицы }

for i:=0 to N-1 do

begin

omega[i]:=CreateSquareMatrix(2);

SetMatrixElement(omega[i],1,1,1);

SetMatrixElement(omega[i],1,2,0);

SetMatrixElement(omega[i],2,1,0);

SetMatrixElement(omega[i],2,2,1);

for j:=0 to i do omega[i]:=MultipleMatrixOnMatrix(matrix_eps[j],omega[i]);

end;

{ решение алгебраической системы двух уравнений относительно коэфф. A и B }

c1.x:=0;

c1.y:=5;

c2.x:=25;

c2.y:=0;

c3.x:=0;

c3.y:=10;

c4.x:=cos(2*k);

c4.y:=sin(2*k);

c5.x:=getmatrixelement(omega[n-1],1,1);

c5.y:=0;

c6.x:=getmatrixelement(omega[n-1],1,2);

c6.y:=0;

c7.x:=getmatrixelement(omega[n-1],2,1);

c7.y:=0;

c8.x:=getmatrixelement(omega[n-1],2,2);

c8.y:=0;

{ подсчет коэффициентов A и B }

Coefficient_A:=cDiv(cSum(cSub(cSum(c7,cMul(c1,c8)),cMul(c1,c5)),cMul(c2,c6)),

cSum(cSum(cSub(cMul(c1,c8),c7),cMul(c1,c5)),cMul(c2,c6)));

Coefficient_B:=cDiv(cMul(c3,cSub(cMul(c5,c8),cMul(c6,c7))),

cMul(c4,cSum(cSub(cMul(c1,c8),c7),cSum(cMul(c1,c5),cMul(c2,c6)))));

{ вывод коэффициентов A и B }

Edit_Re_A.Text:=FloatToStr(Coefficient_A.x);

Edit_Im_A.Text:=FloatToStr(Coefficient_A.y);

Edit_Re_B.Text:=FloatToStr(Coefficient_B.x);

Edit_Im_B.Text:=FloatToStr(Coefficient_B.y);

{ проверка энергетического критерия |A|^2+|B|^2 = 1 }

Edit_Check.Text:=FloatToStr((cNorm(Coefficient_A)+cNorm(Coefficient_B)));

setlength(u,N);

{ вычисление макс и мин функции Re u(x)}

for i:=0 to N-1 do u[i]:=res_u(i*dx,Coefficient_A);

Max_y1:=-1000000; Min_y1:= 1000000;

for i:=0 to N-1 do

begin

if u[i].x> Max_y1 then Max_y1:=u[i].x;

if u[i].x< Min_y1 then Min_y1:=u[i].x;

end;

setlength(u,N);

{ вычисление макс и мин функции Im u(x)}

for i:=0 to N-1 do u[i]:=res_u(i*dx,Coefficient_A);

Max_y2:=-1000000; Min_y2:= 1000000;

for i:=0 to N-1 do

begin

if u[i].y> Max_y2 then Max_y2:=u[i].y;

if u[i].y< Min_y2 then Min_y2:=u[i].y;

end;

{ вывод макc и мин функций}

Edit2.Text:=FloatToStr(Max_y1);

Edit3.Text:=FloatToStr(Max_y2);

Edit4.Text:=FloatToStr(Min_y1);

Edit5.Text:=FloatToStr(Min_y2);

{ построение функции Re u(x)}

dy1:= (PaintBox1.Height-10)/(Max_y1-Min_y1);

dx1:= PaintBox1.Width/(N);

PaintBox1.Canvas.Pen.Color:=clwhite;

PaintBox1.Canvas.Rectangle(0,0,clientwidth,clientheight);

PaintBox1.Canvas.Pen.Color:=clRed;

PaintBox1.Canvas.Pen.Width:=3;

PaintBox1.Canvas.MoveTo(0,5+round(( Max_y1-u[0].x)*dy1));

for i:=1 to N-1 do PaintBox1.Canvas.LineTo(round(i*dx1), 5+round((Max_y1-u[i].x)*dy1));

{ построение функции Im u(x)}

dy2:= (PaintBox1.Height-10)/(Max_y2-Min_y2);

dx2:= PaintBox1.Width/(N);

PaintBox2.Canvas.Pen.Color:=clwhite;

PaintBox2.Canvas.Rectangle(0,0,clientwidth,clientheight);

PaintBox2.Canvas.Pen.Color:=clBlue;

PaintBox2.Canvas.Pen.Width:=3;

PaintBox2.Canvas.MoveTo(0,5+round((Max_y2-u[0].y)*dy2));

for i:=1 to N-1 do PaintBox2.Canvas.LineTo(round(i*dx2),5+round(( Max_y2-u[i].y)*dy2));

end;

end.

10. Результаты численного приближенного метода

Для 50 разбиений:

A = - 0,0611195464652435 + 0,0104851082164756i

B = 0,991383258147422 - 0,11519962826488i

|A|² + |B|² = 0,999957255341792

Для 200 разбиений:

A = - 0,0616466779168554 + 0,0107331622023326i

B = 0,991401340163465 - 0,114925310997946i

|A|² + |B|² = 0,999999958054935

Для 1000 разбиений:

A = - 0,0616801839547689 + 0,0107501248563918i

B = 0,991400263648656 - 0,114915216341547i

|A|² + |B|² = 0,999999999986572

11. Графическое сравнение

Действительная часть Re u(x):

50 разбиений

200 разбиений

1000 разбиений

Мнимая часть Im u(x):

50 разбиений

200 разбиений

1000 разбиений

II. Метод последовательных приближений

12. Теоретическая часть

Известно, что решение уравнения

+ k²u = f(x)

(18)

на всей числовой прямой -   x   удовлетворяющее условию отсутствия волн, приходящих из бесконечности, имеет вид

u(x) = ,

(19)

где G(x, x₁) - функция Грина: решение той же задачи с правой частью f(x) =(x – x₁).

Для уравнения (18) с постоянным коэффициентом k² функция Грина выписывается в явном виде

G(x, x₁) = .

(20)

Следовательно, решение уравнения (18) выписывается в явном виде:

u(x) = .

(21)

Перепишем уравнение (1) в эквивалентной форме

+ k²u = (k² - n(x))u(x).

(22)

Тогда, согласно (21), решение уравнения (22) с правой частью (k² - n(x))u(x) можно выписать в виде

u(x) = . (23)

Так как (k² - n(x))  0 при x  0, и при x  a , то несобственный интеграл в (23) заменяется интегралом в конечных пределах

u(x) = .

(24)

Раскрывая модуль в показателе экспоненты для x  0 и для x  a, получим, что

u(x) = для x  0,

(25)

и

u(x) = для x  a.

(26)

Следовательно, решение (24) содержит лишь волны, уходящие в -  и в  (если зависимость решения от времени выбрана в виде ).

Но исходная постановка задачи содержит также единственную волну , приходящую из -  .

Если взять сумму полей (24) и , то получим представление для полного поля во всей области -   x  

u(x) = + . (27)

В сокращенных обозначениях это уравнение записывается в виде

u = Au + f.

(28)

Здесь А – интегральный оператор , действующий на функцию u(x).

Нетрудно убедиться, что функция u(x), представленная в виде (27), удовлетворяет всем условиям исходной задачи. Действительно, правая часть (27) удовлетворяет уравнению

+ k²u =0

вне слоя [0, a] , и уравнению

+ n(x)u = 0

внутри этого слоя, благодаря свойствам функции Грина (20) .

13. Реализация метода

Краевые условия (2) также выполняются, что проверяется непосредственно. Кроме того при x  0 решение, в силу представления (25), имеет вид

u(x) = + ,

откуда следует, что

А = .

(29)