VPF::Аппроксимация, метод гира. - Форум программистов

shaakvu

Дата 13.4.2009, 13:22 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 1
Регистрация: 13.4.2009

Репутация: нет
Всего: нет

Нахожу методом Рунге-Кутта первые члены уравнений:

//Рунге-Кутт

...

Код


function f1(x,y: real): real;
begin     f1:= b1*e-b1*x-b2*y end;

function f2(x,y: real): real;
begin     f2:= a1*x-a2*y   end;

function rk1(x, y, h: real): real;
var m: array [0..3] of real;
    h_2: real;
begin
     h_2 := 0.5 * h;
     m[0] := f1(x,y);
     m[1] := f1(x + h_2, y + h_2*m[0]);
     m[2] := f1(x + h_2, y + h_2*m[1]);
     m[3] := f1(x + h, y + h*m[2]);
     rk1 := y + h*(m[0] + 2*m[1] + 2*m[2] + m[3])/6;
end;

function rk2(x, y, h: real): real;
var m: array [0..3] of real;
    h_2: real;
begin
     h_2 := 0.5 * h;
     m[0] := f2(x,y);
     m[1] := f2(x + h_2, y + h_2*m[0]);
     m[2] := f2(x + h_2, y + h_2*m[1]);
     m[3] := f2(x + h, y + h*m[2]);
     rk2 := y + h*(m[0] + 2*m[1] + 2*m[2] + m[3])/6;
end;

...

    y1:=0; t:=0.01; y2:=0; x:=0;

     while x <= 0.1 do begin
           y1:= rk1(x, y1, t);
           y2:= rk2(x, y2, t);
           x := x + t;
           end;

           uw[2]:=y1;
           iw[2]:=y2;

...

Теперь, я могу составить и решать эти уравнения методом Гира (второго порядка, по условию задачи):

Код



procedure Raschet(i,R1,R2,C1,L,E,H: extended; Kr,v: integer; var mu1,mi,iw,uw:mass; var Mdl,Mh: MassK);
var G,P,X,Z,Ux,Ix,RaznI,RaznU1,RaznU2,un1,un2,in1,Unm1,Snm1,Snm2,Inm2,Inm1, unew1,inew: extended;
    i2,Otm: integer;

       begin
       uw[i+1]:=(3/4)*uw[i]-(1/3)*uw[i-1]-(2/3)*h*(b1*E-b1*uw[i]-b2*iw[i]);
       iw[i+1]:=(3/4)*iw[i]-(1/3)*iw[i-1]-(2/3)*h*(a1*uw[i]-a2*iw[i]);
       end;

Уравнения составлены правильно, в этом я не сомневаюсь (надеюсь не зря). Далее я вызываю функцию расчёт чтобы посчитать всё:

Код


           for i:= 2 to 100 do // 100 - ограничение на предельное кол-во итераций
           begin
             repeat
             Raschet(i,R1,R2,C1,L,E,H,Kr,v,mu1,mi,iw,uw,Mdl,Mh);
             h:=h/2; // уменьшение шага 
             until abs(uw[i]-uw[i-1])<eps; // пока не будет достигнута искомая погрешность
           end;

И вот этот кусок кода увы не работает, 100 знаков после запятой не предел, хотя eps = 0,001

:'(

Терзаю гугл уже часов 19 на русскоязычных ресурсах про нахождение погрешности ничего, а англоязычный гугл парсит меня на code gear вместо Гира :(

Буду благодарен любой помощи или ссылке - можно на англоязычный ресурс. Прошу извинить за не очень чистый код - переписываю его всю ночь.

M
THandle

Модератор: используйте подсветку кода!

1 Пользователей читают эту тему (1 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема \| Delphi: Общие вопросы \| Следующая тема »