Ирина Эланс

Автор который поможет с любыми образовательными и учебными заданиями

Уравнения в частных производных

Содержание

Постановка задачи. Граничные и начальные условия.…………………..………… 2
Физический смысл уравнения……………………………………………………….. 2
Решение задачи методом Фурье……………………………………………………... 4
Результаты численных экспериментов………………………………………….…… 9
Листинг программы…………………………………………………………………... 17

Постановка задачи.

Методом Фурье решить начальную краевую задачу для волнового уравнения,

в одномерном случае. С условиями: на левом конце – условие Дирихле, на правом конце – условие Неймана.

(1)

-начальные условия (2)

- граничные условия (3)

Физический смысл уравнения.

Уравнение (1) , это волновое уравнение – уравнение колебаний струны.

В одномерном случаи:

При d=1, то уравнение (1) имеет вид : – частный случай. _где- коэффициент упругости

Пусть имеется упругая струна, совершающая колебания под действием внешней силы F.

Функция описывает положение струны в момент времени в точке .

Если d=2, , то уравнение (1) имеет вид :

- уравнение колебания мембраны.

Если d=3, то волновое уравнение описывает:

1) процессы распространения звука в однородной среде;

2) процессы распространения электромагнитных волн в непроводящей среде.

При граничных условиях (3), уравнение (1) описывает колебания струны. У которой один конец закреплен, другой свободен, т.е в любой момент времени t. На свободном конце x=l, натяжение (E- модуль упругости материала струны; S-площадь поперечного сечения струны) равно нулю (нет внешних сил) и, следовательно, .

Решение задачи методом Фурье.

Требуется найти решение уравнения

(1)

удовлетворяющее:

начальным условиям:

(2)

и граничным условиям:

(3)

Будем искать (не равное тождественно нулю) частное решения уравнения (1), удовлетворяющее начальным условиям (2) и граничным условиям (3), в виде

u (x, t) = X (x) T (t).

Подставляя в уравнение (1), получаем:

X (x) T′′(t) = a² X′′(x) T(t).

Разделив переменные получим:

В левой части этого равенства стоит функция, которая не зависит от x, справа – функция, не зависящая от t.

=> левая и правая части не зависят ни от x, ни от t, т. е. равны постоянному числу. Обозначим его через – λ, где λ > 0. Получим следующее уравнение:

Из этих равенств получаем два уравнения:

(3)

(4)

Подставив в уравнение u (x, t) = X (x) T (t) граничные условия, получаем

Откуда получим спектральную задачу Штурма -Лиувиля.

(3)

(4)

(5)

Решаем задачу (3):

Пусть , тогда

т.к

получаем:

Решаем задачу (4):

При , ,

где произвольные постоянные.

Функция имеет вид:

удовлетворяет (1)-(3), значит следующий ряд является решением уравнения:

т. к , то

Это разложение функции в ряд Фурье по:

т.к , то

Это разложение функции в ряд Фурье по:

Осталось подобрать коэффициенты , :

Для нахождения частных решений задачи (1)-(3) остаётся подобрать функции ,.

В качестве , возьмем функций:

тогда коэффициенты примут вид:

или

тогда коэффициенты примут вид:

Невязка(относительная погрешность).

Для проверки нашего приближенного решения используем невязку:

Для вычисления интегралов воспользуемся формулой средних прямоугольников:

и модернизированной формулой средних прямоугольников:

где n-число разбиения отрезка (0,l), h-шаг.

Результаты численных экспирементов.

Тест I.

1)Вычисления интеграла методом средних прямоугольников.

При

h = 0.2

h=0,01

h=0,001

1)Вычисления интеграла модифицированным методом средних прямоугольников.

При

h=0.2

h=0.01

h=0.001

Тест II.

1)Вычисления интеграла методом средних прямоугольников.

При

h = 0.2

h=0.01

h=0.001

1)Вычисления интеграла модифицированным методом средних прямоугольников.

При

h = 0.2

h = 0.01

h = 0.001

Листинг программы.

//----------------------------------------------------------------------------------------------------------- uses

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

Dialogs, StdCtrls, ExtCtrls, TeEngine, Series, TeeProcs, Chart, DbChart, Math;

type

TForm1 = class(TForm)

ComboBox1: TComboBox;

Button1: TButton;

Label1: TLabel;

RadioGroup1: TRadioGroup;

Edit1: TEdit;

Edit2: TEdit;

Edit3: TEdit;

Edit4: TEdit;

RadioGroup2: TRadioGroup;

Memo1: TMemo;

Memo2: TMemo;

Edit5: TEdit;

Edit6: TEdit;

Chart1: TChart;

Series1: TFastLineSeries;

Series2: TFastLineSeries;

procedure FormCreate(Sender: TObject);

function phi(t: real): real;

function psi(t: real): real;

function v(t: real): real;

function vt(t: real): real;

function r1phi: real;

function r1psi: real;

function r2phi: real;

function r2psi: real;

function norm1phi: real;

function norm2phi: real;

function norm1psi: real;

function norm2psi: real;

function residual1: real;

function residual2: real;

function a(k: integer): real;

function b(k: integer): real;

procedure Button1Click(Sender: TObject);

private

{ Private declarations }

public

{ Public declarations }

end;

var

Form1: TForm1;

N: integer;

m: integer;

x: array of real;

alfa: array of real;

beta: array of real;

implementation

{$R *.dfm}

procedure TForm1.FormCreate(Sender: TObject);

begin

ComboBox1.Items.Add('~0.2');

ComboBox1.Items.Add('~0.01');

ComboBox1.Items.Add('~0.001');

ComboBox1.ItemIndex:=0;

RadioGroup1.Items.Add('phi(x)=exp(x)-1-exp(Pi/2)x; psi(x)=x*x');

RadioGroup1.Items.Add('phi(x)=ln(x+1)-2/(Pi+2)x; psi(x)=exp(x)-1');

RadioGroup1.ItemIndex:=0;

RadioGroup2.Items.Add('Метод средних');

RadioGroup2.Items.Add('Модифицированный метод средних');

RadioGroup2.ItemIndex:=0;

end;

// граничные условия

function TForm1.phi(t: real): real;

begin

case RadioGroup1.ItemIndex of

0: result:=exp(t)-1-exp(Pi/2)*t;

1: result:=ln(t+1)-2/(Pi+2)*t;

end;

function TForm1.psi(t: real): real;

begin

case RadioGroup1.ItemIndex of

0: result:=t*t;

1: result:=exp(t)-1;

end;

// приближенное решение и ее производная по времени

function TForm1.v(t: real): real;

var i: integer;

begin

result:=alfa[0]*sin(t);

if m>1 then

for i:=1 to m-1 do result:=result+alfa[i]*sin((2*i+1)*t);

end;

function TForm1.vt(t: real): real;

var i: integer;

begin

result:=beta[0]*sin(t);

if m>1 then

for i:=1 to m-1 do result:=result+(2*i+1)*beta[i]*sin((2*i+1)*t);

end;

//норма в L2(a;b) функции phi(x)

function TForm1.norm1phi: real;

var i: integer;

begin

result:=phi(x[0])*phi(x[0]);

for i:= 1 to N-2 do result:=result+phi(x[i])*phi(x[i]);

result:=sqrt(result);

end;

//норма в С функции phi(x)

function TForm1.norm2phi: real;

var i: integer;

begin

result:=abs(phi(x[0]));

for i:=1 to N-2 do

if abs(phi(x[i]))>result then result:=abs(phi(x[i]));

end;

//норма в L2(a;b) функции psi(x)

function TForm1.norm1psi: real;

var i: integer;

begin

result:=psi(x[0])*psi(x[0]);

for i:= 1 to N-2 do result:=result+psi(x[i])*psi(x[i]);

result:=sqrt(result);

end;

//норма в С функции psi(x)

function TForm1.norm2psi: real;

var i: integer;

begin

result:=abs(psi(x[0]));

for i:=1 to N-2 do

if abs(psi(x[i]))>result then result:=abs(psi(x[i]));

end;

//-------------------------------------------------------------

function TForm1.r1phi: real;

var i: integer;

begin

result:=Power(phi(x[0])-v(x[0]), 2);

for i:=1 to N-2 do result:=result+Power(phi(x[i])-v(x[i]), 2);

result:=sqrt(result);

end;

function TForm1.r1psi: real;

var i: integer;

begin

result:=Power(psi(x[0])-vt(x[0]), 2);

for i:=1 to N-2 do result:=result+Power(psi(x[i])-vt(x[i]), 2);

result:=sqrt(result);

end;

//-----------------------------------------------------------------

function TForm1.r2phi: real;

var i: integer;

begin

result:=abs(phi(x[0])-v(x[0]));

for i:=1 to N-2 do

if abs(phi(x[i])-v(x[i]))>result

then result:=abs(phi(x[i])-v(x[i]));

end;

function TForm1.r2psi: real;

var i: integer;

begin

result:=abs(psi(x[0])-vt(x[0]));

for i:=1 to N-2 do

if abs(psi(x[i])-vt(x[i]))>result

then result:=abs(psi(x[i])-vt(x[i]));

end;

//невязка 1 норма в L

function TForm1.residual1: real;

begin

result:=(r1phi/norm1phi+r1psi/norm1psi)*100;

end;

//невязка 2 (норма в C[a;b])

function TForm1.residual2: real;

begin

result:=(r2phi/norm2phi+r2psi/norm2psi)*100;

end;

//вычисление коэф-ов

function TForm1.a(k: integer): real;

var i: integer;

begin

case RadioGroup2.ItemIndex of

0: begin

result:=0;

for i:=0 to N-2 do

result:=result+phi(x[i])*sin((2*k+1)*x[i])*Pi/2/(N-1);

result:=result*4/Pi;

end;

1: begin

result:=0;

for i:=0 to N-2 do

result:=result+phi(x[i])*(cos((2*k+1)*(i+1)*Pi/2/(N-1))-cos((2*k+1)*i*Pi/2/(N-1)));

result:=result*(-4)/(2*k+1)/Pi;

end;

function TForm1.b(k: integer): real;

var i: integer;

begin

case RadioGroup2.ItemIndex of

0: begin

result:=0;

for i:=0 to N-2 do

result:=result+psi(x[i])*sin((2*k+1)*x[i])*Pi/2/(N-1);

result:=result*4/(2*k+1)/Pi;

end;

1: begin

result:=0;

for i:=0 to N-2 do

result:=result+psi(x[i])*(cos((2*k+1)*(i+1)*Pi/2/(N-1))-cos((2*k+1)*i*Pi/2/(N-1)));

result:=result*(-4)/(2*k+1)/(2*k+1)/Pi;

end;

procedure TForm1.Button1Click(Sender: TObject);

var i: integer;

begin

// количество узлов разбиения промежутка в зависимости от шага

case ComboBox1.ItemIndex of

0: N:=8;

1: N:=157;

2: N:=1571;

end;

//заполняем массив x[i]

SetLength(x, N);

for i:=0 to N-1 do x[i]:=Pi/2/(N-1)*i;

for i:=0 to N-2 do x[i]:=(x[i+1]+x[i])/2;

//заполняем массив alfa[i] и beta[i]

SetLength(alfa, 101);

for i:=0 to 100 do alfa[i]:=a(i);

SetLength(beta, 101);

for i:=0 to 100 do beta[i]:=b(i);

// на график

Chart1.Series[0].Clear;

Chart1.Series[1].Clear;

for i:=1 to 100 do begin //кол-во слагаемых

m:=i;

Chart1.Series[0].Add(residual1, IntToStr(m));

Chart1.Series[1].Add(residual2, IntToStr(m));

end;

{m:=70;

for i:=0 to 157 do begin

Chart1.Series[0].Add(psi(i/100));

Chart1.Series[1].Add(vt(i/100));

end;}

for i:=0 to 10 do Memo1.Lines.Add('');

for i:=0 to 10 do Memo1.Lines[i]:='a'+IntToStr(i)+'='+FloatToStr(RoundTo(a(i), -5));

for i:=0 to 10 do Memo2.Lines.Add('');

for i:=0 to 10 do Memo2.Lines[i]:='b'+IntToStr(i)+'='+FloatToStr(RoundTo(b(i), -5));

Edit1.Text:=FloatToStr(norm1phi);

Edit2.Text:=FloatToStr(norm1psi);

Edit3.Text:=FloatToStr(norm2phi);

Edit4.Text:=FloatToStr(norm2psi);

Edit5.Text:=FloatToStr(residual1);

Edit6.Text:=FloatToStr(residual2);

end;

end.