Процедура решения задачи Коши по МТРК с автоматическим выбором переменного шага интегрирования

ПРИЛОЖЕНИЕ:

ОТВЕТ(задание 3)

Вспомогательные функции:

function fun1(t,y1,y2:real):real; //дано

begin

fun1:=(y2-y1)*t;

end;

function fun2(t,y1,y2:real):real; //дано

begin

fun2:=(y2+y1)*t;

end;

//функция реализующая формулы (12), все обозначения совпадают

function tform1.yh1(h,t0,y10,y20:real):real;

var

k: array [1..4,1..2] of real;

begin

k[1,1]:=h*fun1(t0,y10,y20);

k[1,2]:=h*fun2(t0,y10,y20);

k[2,1]:=h*fun1(t0+(1/2)*h,y10+(1/2)*k[1,1],y20+(1/2)*k[1,2]);

k[2,2]:=h*fun2(t0+(1/2)*h,y10+(1/2)*k[1,1],y20+(1/2)*k[1,2]);

k[3,1]:=h*fun1(t0+(1/2)*h,y10+(1/2)*k[2,1],y20+(1/2)*k[2,2]);

k[3,2]:=h*fun2(t0+(1/2)*h,y10+(1/2)*k[2,1],y20+(1/2)*k[2,2]);

k[4,1]:=h*fun1(t0+h,y10+k[3,1],y20+k[3,2]);

k[4,2]:=h*fun2(t0+h,y10+k[3,1],y20+k[3,2]);

yh1:=y10+(1/6)*(k[1,1]+2*k[2,1]+2*k[3,1]+k[4,1]);

//yh2:=y20+(1/6)*(k[1,2]+2*k[2,2]+2*k[3,2]+k[4,2]);

end;

//построение полинома Эрмита по 2 точкам и 4-ём значениям в этих точках; см. формулу (13)

function tform1.ermit(t0,t1,ft0,ft1,dft0,dft1,trm:real):real;

begin

ermit:=(sqr(trm-t1)/sqr(t0-t1))*(ft0*(1-2*(trm-t0)/(t0-t1))+dft0*(trm-t0))+

(sqr(trm-t0)/sqr(t1-t0))*(ft1*(1-2*(trm-t1)/(t1-t0))+dft1*(trm-t1));

end;

Основные процедуры:

//процедура решения ЗК по МТРК с автоматическим выбором переменного шага интегрирования по заданной величине Е ЛП, поиск значений в точках выдачи результата

procedure TForm1.Button5Click(Sender: TObject);

var

kol,i,j,k,m:integer;

levgran, pravgran,l,t0:real;

y,y2:array [0..1000, 1..2] of real;

sum,h,x,tr:array [0..1000] of real;

begin

kol:=strtoint(edit3.text); //kol-количество точек разбиения

levgran:=strtofloat(edit1.text); //levgran - t0

pravgran:=strtofloat(edit2.text); //pravgran - t0+T

l:=(pravgran-levgran); //длина промежутка = Т

m:=strtoint(edit13.text); //m-кол-во точек выдачи результата

for i:=1 to 1000 do h[i]:=l/kol; //h[i] - начальный шаг интегрирования

stringgrid1.RowCount:=kol+1; //размер таблицы выдачи результата

stringgrid1.Cols[0][0]:=edit4.text; //заполняем первую строку таблицы (это Н.У.З.К.)

stringgrid1.Cols[1][0]:=edit5.text; //заполняем первую строку таблицы (это Н.У.З.К.)

stringgrid1.Cols[2][0]:=edit7.text; //заполняем первую строку таблицы (это Н.У.З.К.)

t0:=strtofloat(edit4.text); //t0

y[0,1]:=strtofloat(edit5.text); //Н.У.

y[0,2]:=strtofloat(edit7.text); //Н.У.

y2[0,1]:=y[0,1]; //Н.У.

y2[0,2]:=y[0,2]; //Н.У.

for i:=1 to m do tr[i]:=t0+(l/m)*i; //фиксируем точки выдачи результатов(tr[i])

i:=0;

x[0]:=t0;

repeat //начало цикла реализующего поиск приближенного решения

i:=i+1;//переход к следующей точке в которой будет происходить поиск приб. реш.

//далее цикл из двух шагов: 1-ый нахождение y[i] ч/з y[i-1] c заданным наперёд начальным шагом инт-я h; 2-ой нахождение y[i] ч/з y[i-1], но с шагом, который обеспечивает точность Е ЛП

for j:=0 to 1 do begin

y[i,1]:=yh1(h[i],x[i-1],y[i-1,1],y[i-1,2]); //находим по МТРК y[i] (П5); см. что такое yh1

y[i,2]:=yh2(h[i],x[i-1],y[i-1,1],y[i-1,2]); //находим по МТРК y[i] (П5); см. что такое yh2

if j<>1 then begin

y2[i,1]:=yh1(h[i]/2,x[i-1],y[i-1,1],y[i-1,2]); //находим по МТРК y2[i] в точке x[i-1]+h/2 (П6)

y2[i,2]:=yh2(h[i]/2,x[i-1],y[i-1,1],y[i-1,2]);

y2[i,1]:=yh1(h[i]/2,x[i-1]+h[i]/2,y2[i,1],y2[i,2]); //находим по МТРК y2[i] в точке x[i-1]+h/2+h/2 (П6)

y2[i,2]:=yh2(h[i]/2,x[i-1]+h[i]/2,y2[i,1],y2[i,2]);

sum[i]:=sqrt(sqr(y2[i,1]-y[i,1])+sqr(y2[i,2]-y[i,2])); //модуль разности для формулы (11)

h[i]:=(h[i]/2)*exp(ln((15*exp(ln(10)*strtofloat(edit12.text))/sum[i]))*(1/4));  //формула (11)

end;

end;

x[i]:=x[i-1]+h[i]; //новая точка в которой нам предстоит искать приближ. значение реш-я

for k:=1 to m do begin //цикл для нахождения значения полинома Эрмита в точке выдачи результата

if (x[i-1]<tr[k])and(tr[k]<x[i]) then begin //проверка, если точка выдачи оказывается м/ду точками разбиения

stringgrid2.Cols[0][k-1]:=floattostr(tr[k]); //заполнение результата в таблицу

stringgrid2.Cols[1][k-1]:=floattostr(ermit(x[i-1],x[i],y[i-1,1], y[i,1],fun1(x[i-1],y[i-1,1],y[i-1,2]),fun1(x[i],y[i,1],y[i,2]),tr[k])); //заполнение результата в таблицу; см. function ermit

stringgrid2.Cols[2][k-1]:=floattostr(ermit(x[i-1],x[i],y[i-1,2], y[i,2],fun2(x[i-1],y[i-1,1],y[i-1,2]),fun2(x[i],y[i,1],y[i,2]),tr[k]));

end;

if tr[k]=x[i] then begin //проверка, если точка выдачи совпадлает с точкой разбиения

stringgrid2.Cols[1][k-1]:=floattostr(y[i,1]);

stringgrid2.Cols[2][k-1]:=floattostr(y[i,2]);

end

end;

until x[i]>t0+l; //завершения основного цикла, происходит в случае выхода за пределы заданного промежутка интегрирования

//Далее: на всякий случай заполним таблицу со значениями во всех точках дробления

stringgrid1.RowCount:=i+1;

for j:=1 to i do

begin

stringgrid1.Cols[0][j]:=floattostr(x[j]);

stringgrid1.Cols[1][j]:=floattostr(y[j,1]);

stringgrid1.Cols[2][j]:=floattostr(y[j,2]);

end;

end;

Рассмотрим процедуру, реализующую нахождение глобальной погрешности.

procedure TForm1.Button1Click(Sender: TObject);

var

kol,i,j:integer;

levgran, pravgran,h,l,t0,max:real;

y,y2,gp,gp2:array [0..10000, 1..2] of real;

sum:array [0..10000] of real;

begin

kol:=strtoint(edit3.text);

levgran:=strtofloat(edit1.text);

pravgran:=strtofloat(edit2.text);

l:=(pravgran-levgran);

h:=l/kol;

stringgrid1.RowCount:=kol+1;

stringgrid1.Cols[0][0]:=edit4.text;

stringgrid1.Cols[1][0]:=edit5.text;

stringgrid1.Cols[2][0]:=edit7.text;

t0:=strtofloat(edit4.text);

y[0,1]:=strtofloat(edit5.text);

y[0,2]:=strtofloat(edit7.text);

y2[0,1]:=y[0,1];

y2[0,2]:=y[0,2];

//всё что выше повторяет предыдущую процедуру

max:=1;

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

//while max>exp(ln(10)*strtofloat(edit11.text)) do begin

//алгоритм нахождения ГП полностью соответствует (П4)

for i:=1 to kol do  //1-ое: ищем значения y в точках xn с шагом h

begin

y[i,1]:=yh1(h,t0+h*(i-1),y[i-1,1],y[i-1,2]);

y[i,2]:=yh2(h,t0+h*(i-1),y[i-1,1],y[i-1,2]);

end;

for i:=1 to kol*2 do //2-ое: ищем значения y в точках xn с шагом h/2

begin

y2[i,1]:=yh1((h/2),t0+(h/2)*(i-1),y2[i-1,1],y2[i-1,2]);

y2[i,2]:=yh2((h/2),t0+(h/2)*(i-1),y2[i-1,1],y2[i-1,2]);

end;

for i:=1 to kol do  //подсчитываем ГП в каждой точке дробления

begin

gp[i,1]:=(y2[i*2,1]-y[i,1])/(1-1/16);

gp[i,2]:=(y2[i*2,2]-y[i,2])/(1-1/16);

gp2[i,1]:=(y2[i*2,1]-y[i,1])/(16-1);

gp2[i,2]:=(y2[i*2,2]-y[i,2])/(16-1);

sum[i]:=sqrt(sqr(gp2[i,1])+sqr(gp2[i,2])); //фиксируем модуль ГП в каждой точке

end;

max:=sum[1];

for i:=1 to kol-1 do //выбираем максимальную ГП

begin

if max<sum[i+1] then begin j:=i; max:=sum[i+1]; end;

end;

stringgrid1.RowCount:=kol+1;

h:=h/2;

kol:=kol*2;

//end;

//далее аналогично предыдущей процедуре

edit6.Text:=floattostr(h*2);

edit8.Text:=floattostr(max);

edit9.Text:=floattostr(t0+h*2*j);

for i:=1 to trunc(kol/2) do

begin

stringgrid1.Cols[0][i]:=floattostr(t0+h*2*i);

stringgrid1.Cols[1][i]:=floattostr(y2[i*2,1]+gp2[i,1]);

stringgrid1.Cols[2][i]:=floattostr(y2[i*2,2]+gp2[i,2]);

end;

end;