|
Nickolay.info. Обучение. Пишем парсер арифметических выражений на Паскале и применяем его для построения интерполяционного сплайна |
Написать парсер - задача несложная. Гораздо сложней написать хороший парсер, но мы пока такой цели и не ставим. В нашем учебном примере ограничимся разбором и арфиметической оценкой строки на Паскале, зависящей от единственного аргумента, обозначенного x. Писать парсер мы будем тоже на Паскале, поэтому обойдемся только стандартными средствами этого языка. По-хорошему стоило бы сделать всё на Си через структуры и указатели (см., например, парсер выражений из этого моего архива, 4 Кб).
В чем, собственно состоит задача?
У нас есть выражение, записанное по правилам языка Паскаль, например, такое:
-1+sin(cos(2*x+-1.25e-02/x))+11.-3-+3-exp(x)-pi/4
Мы хотим сделать с ним всего 2 вещи:
Желательно также, чтоб парсер не был слишком капризен к лишней паре скобок, умел понимать числа в любом формате и вылавливал элементарные математические ошибки вроде деления на ноль. Для удобства стоит оформить парсер в виде отдельного модуля, чтоб потом подключать его к своим программам. Поддержка вычислений вида 1+++x - также не проблема, раз Паскаль такое принимает, а 1+-x и -x+1 тем более могут понадобиться.
Будем действовать самым очевидным образом - сначала разделим строку на отдельные лексемы, потом перепишем эти лексемы в массив, проверим правильность цепочки лексем, наконец, после этого будем вычислять выражение по правилам языка, постепенно сокращая массив, пока на выходе не останется единственное число - результат расчета.
Под лексемой я имею в виду "минимальную единицу" выражения - скобку, знак операции, имя функции или число.
В интерфейсной части модуля опишем следующие глобальные данные:
const
nlex=12; {доступные функции}
lexems: array [1..nlex] of string =(
'sin','cos','arctan','sqr','sqrt','exp','ln','abs','trunc','frac','round',
'pi'
);
max=80; {макс.число лексем}
maxlen=10; {макс.длина лексемы}
var lex:array[1..max] of string[maxlen]; {лексемы}
types:array[1..max] of char; {типы лексем}
numbers:array[1..max] of real; {числа для расчета}
clex:integer; {текущее число лексем}
runerror:string; {ошибка времени исполнения}
Все остальные лексемы, такие как скобки и знаки операций, будем держать прямо в программе. В строке консоли всего 80 символов, поэтому максимальное число лексем принимаем с запасом равным этому значению. Лексемам будут соответствовать типы, которых намечается шесть:
При анализе выражения понадобятся также отдельные метки начала (S) и конца (E) выражения - некоторые лексемы меняют правила, оказываясь в начале или конце выражения. Все эти метки будут помещаться в массив types. Поскольку мы работаем изначально со строкой, получаемые нами числа придется возвращать также в строковые массивы что чревато большими потерями точности при использовании стандартной процедуры str. Поэтому все числа, получаемые при "упаковке" массива лексем будем дополнительно копировать в вещественный массив numbers и брать их для расчета именно оттуда.
Формально правильное выражение не обязательно способно вычислиться - при расчете могут всплыть деление на ноль, извлечение корня из отрицательного числа и другие неприятности. Поэтому введем строку runerror, куда будем писать арифметические ошибки вычисления.
Самое первое, что следует сделать со строкой выражения - убрать лишние пробелы и привести все символы к одному регистру. Следующие 2 функции это и делают, на всякий случай lowcase работает и с кириллицей DOS, понимаемой Паскалем.
function lowcase (s:string):string; {нижний регистр}
var l,i,c:integer;
begin
l:=length(s);
for i:=1 to l do begin
c:=ord(s[i]);
if (c>=ord('A')) and (c<=ord('Z')) or
(c>=ord('А')) and (c<=ord('П'))
then inc (c,32)
else if (c>=ord('Р')) and (c<=ord('Я'))
then inc(c,80);
s[i]:=chr(c);
end;
lowcase:=s;
end;
function alltrim (s:string):string; {убрать пробелы}
var p:integer;
begin
repeat
p:=pos(' ',s);
if p>0 then delete (s,p,1);
until p=0;
alltrim:=s;
end;
Основная функция парсера parse проконтролирует, не пуста ли оставшаяся строка, вызовет функцию проверки и разбора выражения check, и наконец, если выражение верно, вызовет функцию eval, выполняющую расчет. Можно вызывать check и eval независимо друг от друга, если нужно только проверить выражение или вычислить без проверки. Выглядеть функция parse будет следующим образом:
function parse (s:string; x:real; var e:string):real;
var r:real; i:integer;
begin
e:='';
s:=alltrim(lowcase(s));
if length(s)<1 then begin
e:='Ошибка: Строка пуста';
parse:=0; exit;
end;
if check(s,e)=true then begin
for i:=1 to clex do numbers[i]:=0;
if eval(x,r)=true then begin
parse:=r;
exit;
end
else e:='Ошибка вычисления: '+runerror;
end
else e:='Ошибка разбора: '+e;
end;
Параметр s передает строку с выражением, x - аргумент, так как все расчеты делаются в вещественных числах, функция вернет значение типа real, параметр e служит для передачи вызывающей программе сообщений об ошибках.
Остаются check и eval. Проверку формальной правильности выражения проще всего разбить на несколько этапов:
На этапе 2 нужно внимательно учесть, что знаки "+" и "-" могут обозначать как операцию сложения или вычитания, так и знак числа или аргумента x.
function check (var s:string;var e:string):boolean;
{главная функция проверки}
label next_i;
var l,k,k0,i,j,r:integer;
n:real;
s0,s1:string;
begin
l:=length(s);
k:=0;
{контроль скобок}
for i:=1 to l do begin
if s[i]='(' then inc (k)
else if s[i]=')' then dec(k);
if k<0 then begin
e:='Скобки: '+copy(s,i,l-i+1);
check:=false; exit;
end;
end;
if k<>0 then begin
e:='Непарные скобки';
check:=false; exit;
end;
{обработать цепочки знаков +,-}
repeat
k:=0;
inc(k,replace_strings ('++','+',s));
inc(k,replace_strings ('--','+',s));
inc(k,replace_strings ('+-','-',s));
inc(k,replace_strings ('-+','-',s));
until k=0;
{начало разбора на лексемы}
l:=length(s);
s0:='';
for i:=1 to l do s0:=concat(s0,' ');
{выбрать числа}
i:=1;
next_i:
while i<=l do begin
for j:=l downto i do begin
s1:=copy (s,i,j-i+1);
val (s1,n,r);
if r=0 then begin
for k:=i to j do begin
s0[k]:='n';
end;
if (s[i]='+') or (s[i]='-') then begin {лишние + и - из чисел}
if (i>1) and (s[i-1]<>'+') and (s[i-1]<>'-') and (s[i-1]<>'(') then
s0[i]:='o';
end;
i:=j+1;
goto next_i;
end;
end;
inc (i);
end;
{расставить скобки, операции}
k:=0;
for i:=1 to l do begin
case s[i] of
'+','-': if (s0[i]<>'n') then s0[i]:='o';
'*','/': s0[i]:='o';
'(',')': s0[i]:=s[i];
end;
end;
{стандартные функции}
k:=-1;
repeat
k:=pos(' ',s0);
if k>0 then begin
j:=k;
while s0[k]=' ' do begin
inc(k);
if k>l then break;
end;
s1:=copy (s,j,k-j);
k0:=0;
for i:=1 to nlex do begin
if lexems[i]=s1 then begin
for r:=j to k-1 do s0[r]:='f';
k0:=1;
break;
end;
end;
if k0=0 then for r:=j to k-1 do s0[r]:='_';
end
until k=0;
{аргумент}
for i:=1 to l do if (s[i]='x') and (s0[i]='_') then s0[i]:='x';
{...}
Этот алгоритм расставит метки типов в строку s0 той же размерности, что исходная строка s, содержащая выражение. Для приведенного выше выражения, например, получится:
-1+sin(cos(2*x+-1.25e-02/x))+11.-3-+3-exp(x)-pi/4 -- строка s
nnofff(fff(noxonnnnnnnnnox))onnnononnofff(x)offon -- строка s0
Для исключения лишних цепочек плюсов и минусов вызывается функция replace_strings, делающая строковые замены:
function replace_strings (sold,snew:string;var s:string):integer;
{заменить sold на snew в строке s; вернет число замен}
var n,p:integer;
begin
n:=0;
repeat
p:=pos(sold,s);
if p>0 then begin
delete (s,p,length(sold));
insert (snew,s,p);
inc(n);
end;
until p=0;
replace_strings:=n;
end;
Однако, выражения типа -x+1 после нее остаются, почему и необходим дополнительный контроль (блок "{лишние + и - из чисел}").
Функция oper, которая будет вычислять выражения, тоже может получить подобные конструкции после сокращения цепочки лексем, так что постараемся в нее тоже не забыть вставить блок для дополнительного контроля знаков "+" и "-".
Все неверные лексемы мы помечали подчеркиванием в s0, осталось выбрать их и сообщить, если есть ошибка (продолжение функции check):
{оставшиеся цепочки неверны}
for i:=1 to l do if (s0[i]='_') then begin
e:='Неизвестная лексема: '+copy(s,i,l-i+1);
check:=false; exit;
end;
{правила лексем}
if (get_lexem_arrays (s,s0,e)=false) then begin
check:=false; exit;
end;
check:=true;
end;
Теперь все лексемы помечены буквами типов.
На втором этапе своей работы функция check проверяет, выполняются ли формальные правила следования лексем друг за другом, и заодно наполняет данными массивы lex, types и numbers. Все эти действия выполняются путем вызова подпрограммы get_lexem_arrays:
function get_lexem_arrays (s,s0:string;var e:string):boolean;
{получение массива лексем}
var type0,type1:char;
l,start1,end1:integer;
lexem0,lexem1:string;
begin
s0:=s0+'E';
type0:='S';
type1:=s0[1];
l:=length(s);
start1:=1;
clex:=0;
lexem0:='';
while start1<=l do begin
end1:=start1;
if (s0[end1]='(') or (s0[end1]=')') then inc(end1)
else begin
while s0[end1]=type1 do begin
if end1>l then break;
inc(end1);
end;
end;
lexem1:=copy(s,start1,end1-start1);
if length(lexem1)>maxlen then begin
e:='Слишком длинная лексема: '+lexem1;
get_lexem_arrays:=false; exit;
end;
if (check_rool(type0,type1,lexem0,lexem1)=false) then begin
e:='Неверная последовательность действий: '+copy(s,start1,l-start1);
get_lexem_arrays:=false; exit;
end;
inc(clex);
lex[clex]:=lexem1;
lexem0:=lexem1;
types[clex]:=type1;
type0:=type1;
type1:=s0[start1+length(lexem1)];
inc(start1,length(lexem1));
{writeln (lex[clex],' ',types[clex]);}
end;
lexem0:=lexem1;
lexem1:='';
if (check_rool(type0,type1,lexem0,lexem1)=false) then begin
e:='Неверно последнее действие';
get_lexem_arrays:=false; exit;
end;
get_lexem_arrays:=true;
end;
В свою очередь, get_lexem_arrays вызывает функцию check_rool, передавая ей типы и строковые значения каждой пары соседних лексем. Чтобы отследить начало и конец цепочки, они помечаются специальными лексемами S и E. В check_rool заложены правила следования двух соседних лексем, поэтому ее логика довольно наглядна.
function check_rool(pred,cur:char;lexem0,lexem1:string):boolean;
{проверка правил следования лексем}
begin
check_rool:=false;
case pred of
'S': begin
if (cur='f') or (cur='n') or (cur='x') or (cur='(') or
(cur='o') and ((lexem1='-') or (lexem1='+')) then check_rool:=true;
end;
'f': begin
if lexem0<>'pi' then begin
if cur='(' then check_rool:=true;
end
else begin
if (cur='o') or (cur=')') then check_rool:=true;
end;
end;
'n': begin
if (cur='o') or (cur=')') or (cur='E') then check_rool:=true;
end;
'o': begin
if (cur='n') or (cur='x') or (cur='(') or (cur='f') then check_rool:=true;
end;
'x': begin
if (cur='o') or (cur=')') or (cur='E') then check_rool:=true;
end;
'(': begin
if (cur='x') or (cur='f') or (cur='n') or (cur='(') or
(cur='o') and ((lexem1='+') or (lexem1='-')) then check_rool:=true;
end;
')': begin
if (cur=')') or (cur='o') or (cur='E') then check_rool:=true;
end;
end;
end;
Итак, если все нормально, функция check вернула истину, и можно вызывать функцию eval для вычисления.
function eval (x:real;var r:real):boolean; {главная функция оценки}
begin runerror:='';
repeat
r:=skobka (1,x);
if runerror<>'' then begin
eval:=false; exit;
end;
until kol_items ('(')=0;
r:=oper(1,x);
if runerror<>'' then begin
eval:=false; exit;
end;
eval:=true;
end;
Пока в выражении есть нераскрытые скобки, eval будет раскрывать их, обращаясь к рекурсивной функции skobka. Последняя функция, в свою очередь, вызовет функции обслуживания стандартных подпрограмм и арифметических выражений. Попутно функция eval смотрит, не возникло ли глобальной ошибки выполнения runerror и, в случае чего, возвращает значение "ложь", показывающее, что при данном x посчитать формулу нельзя. Когда скобок не осталось, достаточно оценить последнее арифметическое выражение, в которое превратилась цепочка расчетов.
Для оценки количества оставшихся скобочных выражений (непарные и неправильные скобки мы исключили ранее), не мудрствуя лукаво, вызывается функция kol_items.
function kol_items (c:char):integer; {кол-во лексем типа=c}
var i,k:integer;
begin
k:=0;
for i:=0 to clex do if types[i]=c then inc(k);
kol_items:=k;
end;
Функция skobka имеет следующий вид:
function skobka (n1:integer;x:real):real; {разбор круглых скобок}
var n2,k,er:integer; r:real;
begin
n2:=n1+1; k:=0;
repeat
if types[n2]='(' then r:=skobka(n2,x);
if types[n2]=')'then begin
if k=0 then break;
dec(k);
end;
inc(n2);
if n2=clex then break;
until false;
{здесь имеем самые внутренние скобки}
if n1>1 then begin
if types[n1-1]='f' then begin
r:=func(lex[n1-1],n1,x);
end
else begin
n2:=find_oper(n1);
if n2<>0 then r:=oper(n1,x);
end;
end
else begin
n2:=find_oper(n1);
if n2<>0 then r:=oper(n1,x);
end;
if (types[n1]='(') and (types[n1+2]=')') then begin
skobka:=number_skob(n1,x);
exit;
end;
end;
Кроме обработки скобок, она вызывает подпрограмму для вычисления стандартных функций func:
function func(f:string;n1:integer;x:real):real; {стандартные функции}
var r,r2:real;
begin
r:=oper(n1,x);
if f='sin' then begin
r2:=sin(r);
end
else if f='cos' then begin
r2:=cos(r);
end
else if f='arctan' then begin
r2:=arctan(r);
end
else if f='sqr' then begin
r2:=sqr(r);
end
else if f='sqrt' then begin
if r<0 then begin
runerror:='корень из отрицательного числа: '+f;
r2:=0;
end
else r2:=sqrt(r);
end
else if f='exp' then begin
if abs(r)>88 then begin
runerror:='переполнение: '+f;
r2:=0;
end
else r2:=exp(r);
end
else if f='ln' then begin
if r<=0 then begin
runerror:='логарифм от числа<=0: '+f;
r2:=0;
end
else r2:=ln(r);
end
else if f='abs' then begin
r2:=abs(r);
end
else if f='trunc' then begin
r2:=trunc(r);
end
else if f='frac' then begin
r2:=frac(r);
end
else if f='round' then begin
r2:=round(r);
end;
pack_lexems (n1-1,n1+2,r2);
func:=r2;
end;
и подпрограмму поиска арифметических выражений find_oper:
function find_oper (n1:integer):integer; {найти арифм.выражение с поз. n1}
var n2,k:integer;
begin
n2:=n1+1; k:=0;
repeat
if types[n2]=')' then break;
if not ( (types[n2]='n') or (types[n2]='x') or (types[n2]='o') or
(types[n2]='f') and (lex[n2]='pi') ) then inc(k);
inc(n2);
until n2=clex;
if k=0 then find_oper:=n2
else find_oper:=0;
end;
Если в самых внутренних скобках удалось найти подходящее выражение, оно вычисляется функцией oper.
function oper (n1:integer;x:real):real; {арифметические операции}
label again;
var i:integer; r:real;
begin
again:
for i:=1 to clex do begin {убираем "лишние" знаки +,- в числа}
if (types[i]='o') and ((lex[i]='+') or (lex[i]='-')) then begin
if (i=1) and ((types[2]='x') or (types[2]='n'))
or (
(i<clex) and
((types[i+1]='x') or (types[i+1]='n')) and
((types[i-1]='(') or (types[i+1]='o'))
) then begin
pack_number_sign (i,x);
end;
end;
end;
if (lex[n1]='(') and (lex[n1+2]=')') then begin {числа в ()}
oper:=number(n1+1,x);
exit;
end
else begin
i:=n1+1;
while (types[i]<>')') and (i<=clex) do begin
if (lex[i]='*') or (lex[i]='/') then begin {* , /}
do_oper (i,x);
goto again; {могли образоваться новые "висячие" знаки +,-}
end
else inc(i);
end;
i:=n1+1;
while (types[i]<>')') and (i<=clex) do begin
if (lex[i]='+') or (lex[i]='-') then begin {+ , -}
do_oper (i,x);
goto again;
end
else inc(i);
end;
if (lex[i-2]='(') and (lex[i]=')') then begin {осталось (x)}
oper:=number(i-1,x);
exit;
end
else if clex=1 then begin {осталось x}
oper:=number(1,x);
exit;
end;
{-x,+x остаться не могут - убрано выше}
end;
end;
Эта функция "доводит" свое выражение до одного значения. "Тонкая" операция исключения лишнего + и -, объявившегося перед числом после сокращения цепочки лексем, возложена на процедуру pack_number_sign:
procedure pack_number_sign (i:integer;x:real); {убрать знак числа в поз.i}
var k:integer; r:real;
begin
k:=0;
if lex[i]='+' then k:=1
else if lex[i]='-' then k:=-1;
if k<>0 then begin
r:=k*number (i+1,x);
pack_lexems (i,i+1,r);
end;
end;
Непосредственный расчет выражения вида A+B делает подпрограмма do_oper:
function do_oper (n:integer;x:real):real; {вычисление арифметики}
var r:real;
begin
if lex[n]='*' then r:=number(n-1,x)*number(n+1,x)
else if lex[n]='/' then begin
r:=number(n+1,x);
if r=0 then begin
runerror:='деление на 0';
r:=0;
end
else r:=number(n-1,x)/r;
end
else if lex[n]='+' then r:=number(n-1,x)+number(n+1,x)
else if lex[n]='-' then r:=number(n-1,x)-number(n+1,x);
pack_lexems (n-1,n+1,r);
do_oper:=r;
end;
После вычислений цепочку лексем необходимо "паковать", исключая уже обработанное:
procedure pack_lexems (n1,n2:integer;x:real); {упаковать массивы лексем}
var i,k:integer; st:string[maxlen];
begin
if (n1>n2) or (n1<1) or (n2>clex) then exit;
k:=n2-n1;
for i:=n2+1 to clex do begin
lex[i-k]:=lex[i];
types[i-k]:=types[i];
numbers[i-k]:=numbers[i];
end;
str(x:maxlen,st);{неточно}
numbers[n1]:=x;{поэтому делаем копию}
lex[n1]:=st;
types[n1]:='n';
dec(clex,k);
end;
Заодно pack_lexems борется с неточным сохранением чисел стандартной процедурой str, дублируя эти числа в массив numbers.
Наконец, функция number_skob нужна для обработки оставшегося в скобках единственного аргумента
function number_skob (n1:integer;x:real):real; {число в скобках}
var r:real;
begin
r:=number(n1+1,x);
pack_lexems (n1,n1+2,r);
number_skob:=r;
end;
а функция number - просто аргумента без скобок:
function number (n1:integer;x:real):real; {обработка числа}
var r:real; er:integer;
begin
if (types[n1]='f') and (lex[n1]='pi') then r:=pi
else if types[n1]='x' then r:=x
else begin
if abs(numbers[n1])>1E-8 then r:=numbers[n1] {точная копия числа}
else begin
val(lex[n1],r,er);
if er<>0 then begin
runerror:='нет числа:'+lex[n1];
r:=0;
end;
end;
end;
number:=r;
end;
При этом функция берет ранее сохраненные в массиве numbers значения.
Остается соединить все в один модуль и получится файл parser.pas, который можно скачать по этой ссылке:
парсер выражений на Паскале (архив *.zip со всеми файлами .pas из статьи и скомпилированным parser.tpu, 12 Кб).
Для тестирования модуля напишем маленький файл parse_me.pas.
uses crt,parser;
var s,e:string;
r,rr,x:real;
begin
clrscr;
s:='-1+sin(cos(2*x+-1.25e-02/x))+11.-3-+3-exp(x)-pi/4';
x:=1;
r:=parse(s,x,e);
if e<>'' then writeln (e)
else writeln ('Result=',r);
rr:=-1+sin(cos(2*x+-1.25e-02/x))+11.-3-+3-exp(x)-pi/4;
writeln ('Control=',rr);
reset (input); readln;
end.
Надеюсь, Вы не забыли, что для подключения модуля его нужно иметь скомпилированным в одной из папок, прописанных в настройке EXE & TPU Directory меню Options.Directories Турбо Паскаля. Скомпилировать исходник модуля на диск можно, если открыть в редакторе файл parser.pas, в меню Compile.Destination выставить значение Disk и выполнить команду меню Compile.Build.
Вместо содержимого строки s и выражения, присвоенного переменной rr, можно поставить и любое другое выражение, например,
-(-pi/4+trunc(x/8)++-1.5*(x-11/x))
тоже прекрасно сработает. Кстати, и
((-(((-pi/4+trunc(x/8)++-1.5*(x-11/x))))))
теперь наш парсер не смутит.
На практике удобнее один раз проверить с помощью check правильность заданной пользователем функции, а потом уже вычислять ее для разных значений x с помощью eval. Правда, тогда придется снабдить наш модуль дополнительной функцией, возвращающей вызывающей программе данные разбора, такие как переменная clex, массивы lex и types - мы-то все делаем в глобальных переменных и массивах. Но эту задачу я возлагаю на Вас, а в примере ниже буду просто предварительно проверять правильность введенной функции с помощью check, а потом вызывать parse, которая повторит всю ее трудоемкую работу. Лучше было бы вызвать eval напрямую.
Чтобы не быть голословными, давайте построим любимый мной с детства интерполяционный сплайн. Это такая гладкая кривая с неперерывной второй производной, проходящая через набор заданных точек.
Далее приводится листинг файла spline.pas, содержащего подпрограммы для построения интерполяционного сплайна, а также подключающего парсер для того, чтобы построить сплайн по таблице значений произвольной функции. В программе не делаются все проверки корректности, их легко добавить. Обратите внимание на вызов парсера из функции f, а также на ввод строки s в главной программе модуля.
{Кубический интерполяционный сплайн}
program Spline;
uses crt,parser; {модуль управления экраном и парсер}
type vector=array [0..100] of real; {Нумеруем точки с нуля}
var x,y,c:vector;
x0,x9,h,x1,p,p1,p2,e:real;
n,i:integer;
s,er_str:string;
function f(s:string;x:real;var er_str:string):real; {функция}
var r:real;
begin
er_str:='';
r:=parse(s,x,er_str);
if er_str='' then f:=r
else begin
writeln (er_str);
f:=0;
end;
end;
procedure MyInput (var n:integer; var x,y:vector); {Ввод исходных данных}
var i:integer;
begin
writeln ('Введите границы по оси X для построения полинома:');
read (x0,x9);
writeln ('Введите шаг по X для построения значений полинома:');
read (h);
n:=round((x9-x0)/h+1);
x[0]:=x0; y[0]:=f(s,x0,er_str);
for i:=1 to n do begin {Вычислить значения функции}
x[i]:=x[i-1]+h;
y[i]:=f(s,x[i],er_str);
end;
end;
procedure Coeff(n:integer; var x,f,c:vector);
{Вычисление коээфициентов сплайна}
var i,j,m:integer;
a,b,r:real;
k:vector;
begin
{Прямой ход прогонки}
k[1]:=0; c[1]:=0;
for i:=2 to n do begin
j:=i-1;
m:=j-1;
a:=x[i]-x[j];
b:=x[j]-x[m];
r:=2*(a+b)-b*c[j];
c[i]:=a/r;
k[i]:=(3.0*((f[i]-f[j])/a-(f[j]-f[m])/b)-b*k[j])/r;
end;
{Обратный ход прогонки}
c[n]:=k[n];
for i:=n-1 downto 2 do c[i]:=k[i]-c[i]*c[i+1];
end;
procedure Spl (n:integer; var x,f,c:vector; x1:real; var p,p1,p2:real);
{Построение сплайна. x,f - исходные данные, c - вектор коэффициентов,
наденный процедурой Coeff, x1 - значение x, для которого строим сплайн,
p - значение сплайна в точке, p1,p2 - 1-я и 2-я производные}
var i,j:integer;
a,b,d,q,r:real;
begin
i:=1;
while (x1>x[i]) and (i<>n) do i:=i+1; {Ищем номер соседнего узла}
{Промежуточные переменные и коэффициенты}
j:=i-1; a:=f[j]; b:=x[j]; q:=x[i]-b;
r:=x1-b; p:=c[i]; d:=c[i+1];
b:=(f[i]-a)/q - (d+2*p)*q/3.0;
d:=(d-p)/q*r;
{Считаем значения сплайна и его производных:}
p1:=b+r*(2*p+d);
p2:=2*(p+d);
p:=a+r*(b+r*(p+d/3.0))
end;
begin
clrscr; {очистить экран}
writeln ('Построение кубического интерполяционного сплайна');
repeat
writeln ('Введите функцию в виде выражения на Паскале, например, sqr(sqr(x))+5*sqr(x)-x+1');
readln (s);
if (IoResult<>0) or (check(s,er_str)=false) then begin
if er_str<>'' then writeln (e);
writeln ('Введенное выражение неверно, пожалуйста, повторите');
end
else break;
until false;
MyInput (n,x,y);
Coeff (n,x,y,c); {Нашли коэффициенты C с помощью прогонки}
{Строим таблицу значений сплайна}
writeln ('Значение X':19,'Функция':19,'Сплайн':19,'Отн.ошибка':19);
for i:=0 to n do begin
Spl (n,x,y,c,x[i],p,p1,p2);
if abs(y[i])>1E-9 then e:=abs(y[i]-p)/abs(y[i])
else e:=0;
writeln (x[i]:19:8,y[i]:19:8,p:19:8,e:19:8);
end;
writeln ('Enter для выхода...');
reset (input); readln;
end.
Вот скриншот тестового прогона программы spline.pas:
Введите функцию в виде выражения на Паскале, например, sqr(sqr(x))+5*sqr(x)-x+1
sqr(sqr(x))+5*sqr(x)-x+1
Введите границы по оси X для построения полинома:
0 2
Введите шаг по X для построения значений полинома:
0.2
Значение X Функция Сплайн Отн.ошибка
0.00000000 1.00000000 1.00000000 0.00000000
0.20000000 1.00160000 1.00160000 0.00000000
0.40000000 1.42560000 1.42560000 0.00000000
0.60000000 2.32960000 2.32960000 0.00000000
0.80000000 3.80960000 3.80960000 0.00000000
1.00000000 6.00000000 6.00000000 0.00000000
1.20000000 9.07360000 9.07360000 0.00000000
1.40000000 13.24160000 13.24160000 0.00000000
1.60000000 18.75360000 18.75360000 0.00000000
1.80000000 25.89760000 25.89760000 0.00000000
2.00000000 35.00000000 35.00000000 0.00000000
2.20000000 46.42560000 46.42560000 0.00000000
Enter для выхода...
Если программа найдет ошибку времени исполнения, соответствующее значение f(x) будет вычислено как 0. Например, такое произойдет в первой точке при вводе функции 1-2/x, интервала от 0 до 1 и шага 0.5.
 гостевая; E-mail
|