Вычислит.мат_лаб.р

procedure kvadratur(fun:string;a,b,ep

s:real;var res,d::real;)

c=sqrt(3 / 5) h1=(b-a) / 2 c1=c*h1 x2=(b+a) / 2

f1=execute(fun,x2-c1) f3=execute(fun,x2+c1) s1=h1*(5*f1+8*execute(fu n,x2)+5*f3) / 9

n=2

N

N=N*2

h=(b-a)/ n h1=a/2 c1=c*h1 x2=a+h1 x1=x2c1 x3=x2+c1 s2=0

I

I>N

s2=s2+5*execute(fun,x1)+

8*execute(fun,x2)+5*exec

ute(fun,x3)

x1=x1+h

x2=x2+h

x3=x3+h

I=I+1

I

s2=s2*h1/ 9 d=abs(s1-s2)/ 63 s1=s2

n=n*2

d<eps N

res=s2

Возврат

Рисунок 34 - Схема алгоритма квадратурной формулы Гаусса

Варианты заданий для решения задач численного интегрирования и дифференцирования приведены в таблице 7.

81

Дифференцирование с помощью сплайнов

Входные параметры: a,b – интервал дифференцирования; x – точка дифференцирования; n – число точек дифференцирования; fun – вид функции.

Выходные параметры: y_1,y_2 – значения первой и второй производной в точке x соответственно.

Схема алгоритма показана на рисунке 35.

Пример. Продифференцировать функцию e x2 , количество точек разбиения 5, значение точки х=1.

Текст программы:

procedure Spl_Integr(fun:string;a,b,x:real;n:integer;var y_1,y_2:real); var h,h1,h2,aa:real;

y:array of real; begin h:=(b-a)/n; SetLength(y,n-1); aa:=a;

for i:=0 to length(y)-1 do begin

y[i]:=execute(fun,aa); form3.Memo1.Lines.Add(floattostr(y[i])); aa:=aa+h;

end; i:=trunc((x-a)/h+h/2);

h1:=2*h; h2:=h*h; if i=0 then begin

y_1:=(-3*y[0]+4*y[1]-y[2])/h1; y_2:=(2*y[0]-5*y[1]+4*y[2]-y[3])/h2; end;

if (i>0)and(i<n) then begin

y_1:=(-y[i-1]+y[i+1])/h1; y_2:=(y[i-1]-2*y[i]+y[i+1])/h2; end; if i=n then

begin

y_1:=(y[n-2]-4*y[n-1]+3*y[n])/h1; y_2:=(-y[n-3]+4*y[n-2]-5*y[n- 1]+2*y[n])/h2;

end;

end;

Вычисления по программе привели к следующим результатам: Первая производная:3,808 Вторая производная:10,00

82

procedure spl_integr(fun:string;a,b,x: real;n:integer;var y_1,y_2::real;)

h=(b-a)/ n

Setlength(y,n-1)

aa=a

I

I>length(y)-1

y[i]=execute(fun,aa)

aa=aa+h

I=I+1

I

i=trunc((x-a) / h+h/ 2) h1=2*h h2=h*h

+

i=0

y_1=(- 3*y[0]+4*y[1]- y[2])/ h1 y_2=(2*y[0]-5*y[1]+4*y[2]- y[3])/ h2

+

i>0 andi<n

y_1=(- y[i-1]+y[i+1])/ h1 y_2=(y[i-1]-2*y[i]+y[i+1])/ h2

i=n

+

y_1=(y[n-2]-4*y[n-1]+ 3y[n])/ h1 y_2=(-y[n-3]+4*y[n-2]- 5*y[n-1]+2*y[n])/ h2

Y_1

Y_2

Возврат

Рисунок 35 - Схема алгоритма дифференцирования с помощью сплайнов

Варианты заданий приведены в таблице 7.

83

Лабораторная работа №11 Жесткие задачи для систем ОДУ

Метод Гира

Входные параметры: n – порядок системы; y0 – массив из n чисел, содержащий решение в начальной точке; eps – интегрирования; masq – массив размерности n n, содержащий матрицу системы.

Выходные параметры: res – массив из n чисел, содержащий решение системы.

Схема алгоритма показана на рисунке 36. Пример. Решить задачу Коши для системы

y1 11y1 9 y2 y2 9 y1 11y2

при начальных условиях у1(0)=1, у2(0)=0.

Текст программы:

procedure gira(a,b:real;n,kolfun:integer;x:real;var y_1:TFunZnach); var fun:array[1..NumFun,1..4] of real;

h:real; begin h:=(b-a)/n;

for i:=1 to 4 do

begin runge_ku(a,b,1,1,x,y_1); x:=x+h; for k:=1 to kolfun do fun[k,i]:=y_1[k];

end;

for i:=5 to n do begin

for k:=1 to kolfun do begin

y_1[k]:=fun[k,4]; runge_ku(a,b,1,1,x,y_1); y_1[k]:=1/25*(48*fun[k,4]-36*fun[k,3]+16*fun[k,2]-

3*fun[k,1]+12*h*f(y_1[k],x));

fun[k,1]:=fun[k,2]; fun[k,2]:=fun[k,3]; fun[k,3]:=fun[k,4]; fun[k,4]:=y_1[k];

end;

x:=x+h;

end;

end;

Вычисления по программе привели к следующим результатам:

84

Х	У(1)	у(2)
.10000Е	.48703Е	.34270Е
+00	+00	+00
.20000Е	.35432Е	.32700Е
+00	+00	+00
.30000Е	.28565Е	.27417Е
+00	+00	+00
.40000Е	.23483Е	.22550Е
+00	+00	+00
.50000Е	.19396Е	.18492Е
+00	+00	+00
.60000Е	.16060Е	.15159Е
+00	+00	+00
.70000Е	.13330Е	.12430Е
+00	+00	+00
.80000Е	.11095Е	.10195Е
+00	+00	+00
.90000Е	.83650Е-	.82750Е-
+00	01	01
.10000Е	.68668Е-	.67768Е-
+01	01	01
.11000Е	.56402Е-	.55502Е-
+01	01	01
.12000Е	.46359Е-	.45459Е-
+01	01	01
.13000Е	.38137Е-	.37237Е-
+01	01	01
.14000Е	.31405Е-	.30505Е-
+01	01	01
.15000Е	.25894Е-	.24994Е-
+01	01	01
.16000Е	.21381Е-	.20481Е-
+01	01	01
.17000Е	.17687Е-	.16787Е-
+01	01	01
.18000Е	.14662Е-	.13762Е-
+01	01	01
.19000Е	.12185Е-	.11285Е-
+01	01	01
.20000Е	.92579Е-	.91678Е-
+01	02	02

85

procedure gira(a,b:real;n,kolfun:integ er;x:real;var y_1:TFunZnach)

h=(b-a)/ n

I

I>4

runge_ku(a,b,1,1,x, y_1)

x=x+h

K

K>kolfun

fun[k,i]=y_1[k]

K=K+1

K

I=I+1

I

I

I>5

K

K>kolfun

y_1[k]=fun[k,4]

runge_ku(a,b,1,1,x, y_1)

y_1[k]=1/25*(48*fun[k,4]-3 6*fun[k,3]+16*fun[k,2]-3*fu n[k,1]+12*h*f(y_1[k],x))

fun[k,1]=fun[k,2]

fun[k,2]=fun[k,3]

fun[k,3]=fun[k,4] fun[k,4]=y_1[k]

K=K+1

K

x=x+h

I=I+1

I

Возврат

.

Рисунок 36 - Схема алгоритма метода Гира

86

Варианты заданий для решения Жестких задач для систем ОДУ приведены в таблице 8.

Метод Ракитского (матричной экспоненты).

Входные параметры: n – порядок системы; x – начальное значение x;

h – шаг интегрирования h; Q – массив из n чисел, содержащий решение в начальной точке; F – массив размерности n n, содержащий матрицу системы.

Выходные параметры: X – массив из n чисел, содержащий решение системы.

Схема алгоритма показана на рисунке 37. Пример. Решить задачу Коши для системы

y1 11y1 9 y2 y2 9 y1 11y2

при начальных условиях у1(0)=1, у2(0)=0.

Текст программы:

Procedure MetRak(N:Word;Var H:Real;A:MasReal;Var

F,Q:MasReal;Nts,Ind:Integer);

Var

RabMas1: MasReal;
RabMas2: MasReal;
S,Fak : Real;
i,j,k : Byte;
BEGIN
IF Ind=0 Then
Begin
S:=0;
For i := 1 To N Do
For j := 1 To N Do	S:=S+A[i,j]*A[i,j];
S:=Sqrt(S); H:=0.1/S;		Nts := 0; RabMas1 := A;
For I := 1 To N Do
Begin
For J := 1 To N Do	Begin
Q[i][j] := A[i][j]*H/2;		F[i][j] := A[i][j]*H
End;
Q[i][i] := Q[i][i] + 1;		F[i][i] := F[i][i] + 1
End;
Fak := 1;
For K := 2 To 4 Do
Begin Fak := Fak*k;		MulMat(N, RabMas1, A, RabMas2);

87

For I := 1 To N Do

For J := 1 To N Do Begin

F[i][j] := F[i][j] + RabMas2[i][j]*Step(H, k)/Fak;

IF k < 4 Then Q[i][j] := Q[i][j] + RabMas2[i][j]*Step(H, k)/(Fak*(k+1)); End;

RabMas1 := RabMas2 End;

MulConst(N, h, Q) End;

IF Nts > 0 Then

Begin

H := H*Step(2,Nts); MulMat(N, F, F, RabMas1);

RabMas2 := F; F := RabMas1;

For I := 1 To N Do RabMas2[i][i] := RabMas2[i][i] + 1;

MulMat(N, Q, RabMas2, RabMas1);

Q := RabMas1

End;

End;

Вычисления по программе привели к следующим результатам:

88

х	у(1)	у(2)
.10000Е	.47703Е	.34170Е
+00	+00	+00
.20000Е	.34432Е	.32600Е
+00	+00	+00
.30000Е	.27565Е	.27317Е
+00	+00	+00
.40000Е	.22483Е	.22450Е
+00	+00	+00
.50000Е	.18396Е	.18392Е
+00	+00	+00
.60000Е	.15060Е	.15059Е
+00	+00	+00
.70000Е	.12330Е	.12330Е
+00	+00	+00
.80000Е	.10095Е	.10095Е
+00	+00	+00
.90000Е	.82650Е-	.82650Е-
+00	01	01
.10000Е	.67668Е-	.67668Е-
+01	01	01
.11000Е	.55402Е-	.55402Е-
+01	01	01
.12000Е	.45359Е-	.45359Е-
+01	01	01
.13000Е	.37137Е-	.37137Е-
+01	01	01
.14000Е	.30405Е-	.30405Е-
+01	01	01
.15000Е	.24894Е-	.24894Е-
+01	01	01
.16000Е	.20381Е-	.20381Е-
+01	01	01
.17000Е	.16687Е-	.16687Е-
+01	01	01
.18000Е	.13662Е-	.13662Е-
+01	01	01
.19000Е	.11185Е-	.11185Е-
+01	01	01
.20000Е	.91579Е-	.91578Е-
+01	02	02

89

90