Библиотека численных методов для PascalABC.NET

Среда PascalABC.NET
#21

Ну, неэффективно!

Вот как надо эффективнее:

  class function operator -(a, b: fraction): fraction;
  begin
    var n := GCD(a.denominator, b.denominator);
    Result := new fraction(a.numerator * (b.denominator div n)- b.numerator * (a.denominator div n),
        a.denominator div n * b.denominator);
  end;      
...
  constructor Create(num: BigInteger := 0; denom: BigInteger := 1);
  begin
    var n := GCD(num,denom);
    numerator := num/n;
    denominator := denom/n
  end;

И убрать все эти внешние Trunc.

И += -= переопределить

#22

? При компиляции рекурсивно десериализуются в дерево только используемые функции. Все крайне эффективно.

1 лайк
#23

Безусловно, так лучше. Сразу в голову не пришло, я же спешил проверить))) += и прочие тоже можно переопределить, это не проблема.

Так что, идем этим путем?

#24

Ну да, идём.

И a.SwapRows(i,k);

И протестировать бы

#25

Я попробовал, но Вы же знаете аксиому: тестирование никогда не следует поручать разработчику кода.

#26

Да ладно, бросьте

#27

Я проверил поведение на 2,3 и 10 точках для полинома 3й степени и 4й. Сплайн их прекрасно аппроксимирует. А быстро осциллирующие функции с разрывами - это за границами теории сплайнов. Там же идея проста, как мычание: по каждым четырем соседним точкам находим коэффициенты единственного существующего полинома третьей степени, что дает гладкую трижды дифференцируемую на этом интервале функцию. И этот полином, как лекало, скользит по точкам, плавно их сшивая. Кстати, зная формулу кубического полинома :grinning:, мы можем его аналитически дифференцировать и во всем заданном диапазоне интерполировать не только значение функции, но и весьма точное значение первых и вторых производных. А это тоже неплохо, потому что “выбрыки” численного дифференцирования уже давно стали притчей во языцех.

#28

Ну и вставьте ваши тесты в виде Assertов в конце программы

#29

А сакральный смысл этого?

#30

Ну, мы поверим, все начнут пользоваться, всемирная известность и всё такое

2 лайка
#31

Так их потом отрезать придется же, в библиотеке класса совершенно нечего делать тестированию. Или Вы о том, чтобы к библиотеке был внешний модуль, реализующий тестирование и на нем можно было проверять работу библиотеки на очередной версии компилятора?[quote=“Admin, post:258, topic:143, full:true”] всемирная известность и всё такое [/quote] Я за свою жизнь уже столько наелся известности, что временами хотелось спрятаться))) Слава богу, это уже в прошлом. Знаете, когда в свое время уже из Минрадиопрома нашему руководителю звонили и говорили, что в Баку в Минлегпроме c EC-1036 проблема, не могут второй месяц решить, просят Вашего специалиста командировать, а я сижу с другой командировкой на руках - такая известность не шибко радует.

#32

Ну да, хотелось бы

#33

Это несложно, сделаю, конечно.

#34

Примерно так?

uses StudLib;

begin
  var n:=10;
  var x:=ArrGen(n,i->real(i+1));
  var y:=ArrGen(n,i->sqr(x[i])*x[i]);
  var Sp:=new Spline(x,y); // создаем сплайн с заданными узлами интерполяции.
  var r:=Sp.Value(2.5);
  Assert(r=15.625,'Ожидалось 15.625, получено '+r);
  r:=Sp.Value(1);
  Assert(r=1,'Ожидалось 1, получено '+r);
  r:=Sp.Value(10);
  Assert(r=1000,'Ожидалось 1000, получено '+r);
  y:=ArrGen(n,i->Power(x[i],5));
  Sp:=new Spline(x,y);
  r:=Sp.Value(2.5);
  var eps:=1e-10; // может оказаться зависимой от OC
  Assert(abs(r-96.2230762150221)<=eps,'Ожидалось 96.2230762150221, получено '+r);
  r:=Sp.Value(1);
  Assert(r=1,'Ожидалось 1, получено '+r);
  r:=Sp.Value(10);
  Assert(abs(r-94936.1745213549)<=eps,'Ожидалось 94936.1745213549, получено '+r);
  y:=ArrGen(n,i->(3*Power(x[i],-4)-7*x[i])*sin(x[i]/0.7));
  Sp:=new Spline(x,y);
  r:=Sp.Value(2.5);
  Assert(abs(r-8.18890114377246)<=eps,'Ожидалось 8.18890114377246, получено '+r);
  r:=Sp.Value(1);
  Assert(abs(r+3.9596123054885)<=eps,'Ожидалось -3.9596123054885, получено '+r);
  r:=Sp.Value(10);
  Assert(abs(r+579.522280926788)<=eps,'Ожидалось -579.522280926788, получено '+r);
end.

А вот это я не понял…

#35

Да, примерно так. Только непонятно, откуда берутся числа.

Мне казалось, что в одном месте Вы меняете циклом строки. Может, я ошибался

Ещё один момент. А почему бы не сделать инициализацию сплайна набором точек, а не двумя массивами?

#36

Числа один раз вычисляются и подставляются в код. Это же интерполяция, тут не существует возможности получить точное значение “с помощью карандаша и бумаги” за исключением случаев, когда исходная функция является полиномом степени не выше 3. Отсюда и “эпсилон”, появившейся при аппроксимации сплайном полинома пятой степени и колебательной функции.

Я посчитал, что обычные массивы начинающим могут быть более понятны, чем массив из кортежа, да и получить их проще. Вспомнил, что в при работе с графикой, например в полилинии, получать массив точек не всегда удобно. И еще, поскольку библиотека будет содержать множество программ, часть из них будут использовать не только координаты точек и были мысли о некотором единообразии данных, посему также победили массивы.

Мы могли бы работать с понятием точки, если бы Вы ввели в ядро такой примитив и дали набор функций для манипулирования точками.

Сделал тест покомпактнее.

uses StudLib;

procedure TestSpline(x:real; F:real->real; S:Spline);
begin
  var r1:=F(x);
  var r2:=S.Value(x);
  Assert(r1=r2,'Ожидалось '+r1+', получено '+r2);
end;

procedure TestSpline(x,y,eps:real; S:Spline);
begin
  var r:=S.Value(x);
  Assert(Abs(r-y)<=eps,'Ожидалось '+y+', получено '+r);
end;

begin
  var n:=10;
  var x:=ArrGen(n,i->real(i+1));
  var y:=ArrGen(n,i->sqr(x[i])*x[i]);
  var Sp:=new Spline(x,y); // создаем сплайн с заданными узлами интерполяции.
  TestSpline(1,x->x*x*x,Sp); // левая точка
  TestSpline(2.5,x->x*x*x,Sp); // внутри
  TestSpline(10,x->x*x*x,Sp); // правая точка
  //
  y:=ArrGen(n,i->Power(x[i],5));
  Sp:=new Spline(x,y);
  TestSpline(1,1,real.Epsilon,Sp);
  TestSpline(2.5,96.2230762150221,real.Epsilon,Sp);
  TestSpline(10,94936.1745213549,1e-10,Sp);
  //
  y:=ArrGen(n,i->(3*Power(x[i],-4)-7*x[i])*sin(x[i]/0.7));
  Sp:=new Spline(x,y);
  TestSpline(1,-3.9596123054885,1e-14,Sp);
  TestSpline(2.5,8.18890114377246,1e-14,Sp);
  TestSpline(10,-579.522280926788,1e-12,Sp);
  //
  Writeln('Проверка завершена')
end.
#37

А Вы их вычисляете на бумаге? Нельзя ли вычисления тоже приводить в тесте?

Да бросьте. Опять-таки, тут подменяются понятия - начинающим программировать или начинающим программировать численные методы. Начинающие программировать - это школьники, нередко 5-7 классов. Начинающие программировать численные методы - это студенты 2-3 курсов.

Почему?

uses GraphABC;

begin
  Polyline(Arr((10,10),(10,20),(20,20)));
  Polyline((10,10),(10,20),(20,20));
end.

Не буду - это примитив предметной области. Введите его у себя.

#38

По-моему, массив из кортежей - это обычный массив

Почему?

  var a := Partition(1.0,2.0,10).Select(x->Point(x,f(x)));
#39

Конечно же нет, я их вычислил для случая полинома пятой степени и синусоидальной функции при помощи этой же программы. Результат сплайн-интерполяции при помощи конкретной программы на бумаге не вычислишь: на то он и численный метод. Мы же не оцениваем результат, как отклонение от истинного значения, а проверяем, работает ли модуль, т.е. эти тесты - просто верификация факта, что в очередной версии PascalABC.NET библиотека “не развалилась”. Там даже если аналитически считать, мы сразу упираемся в: а) конечную точность вычисления значения функции б) конечную точность решения системы уравнений. А вот полином третьей степени - там именно вычисление функции сделано.

Т.е. Вы предлагаете сделать класс Point со своим конструктором, определить над ним набор операций и т.д.? Ведь давайте на минуту представим, что интерполяция - лишь часть задачи и исходные точки (может быть, отдельно их абсциссы даже) нужны также и в других целях. Делать Point.UnZip? И еще. Если график построить захочется, придется подключать GraphABC. А там Point уже есть…

#40

А главное, я пока не увидел явной пользы от этого. Коды реализации численных методов это только усложнит. Практически там нет операторов, где одновременно используются абсцисса и ордината точки, а вместо x[i] и y[i] придется писать что-то в стиле P[i][0] и P[i][1]? При этом резко рушится информативность кода. А если сделать свойства и писать P[i].x, P[i].y, то вроде и кортеж непонятно зачем, а код просто получается длиннее. Можно, конечно, принимать параметры массивом кортежей и тут же рассыпать его в два одномерных массива х и у, но это лишь двойной расход памяти. Может, я чего-то не понял, просветите.