SDev.Pro - разработка на заказ

http://sdev.pro - это:

1) Разработка программного обеспечения

Мы используем платформы ASP.NET MVC, LAMP, Atmel, Android, iOS, DirectShow и базы данных MS SQL, Oracle, PostgreSQL, а также облака на MS Azure и Amazon AWS для реализации любых Ваших идей.

2)Внедрение программного обеспечения

Мы предлагаем настройку и доработку решений на базе систем электронного документооборота MS SharePoint Server, геоинформационных систем на GeoServer, систем отчетности на MS SQL Reporting Services, а также облачных решений на MS Office 365.

http://sdev.pro

Sep 3, 2009

Наука: интерполяция сплайнами

Выкладываю пример (Delphi) интерполяции сплайнами. Основан на примерах из книги Мудров А.Е. “Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль (Томск, 1991)”

// Пример интерполяции сплайнами
// Игорь Подсекин (wonderu@gmail.com)
// Лицензия LGPL
//
// http://dotnet.wonderu.com
// http://directshow.wonderu.com

program spline;

{$APPTYPE CONSOLE}

type
vec = array[0..100] of Real;

var
i, k, n: Integer;
x1, x0, x9, h, p, p1,p2:Real;
x, f, c: vec;

//Формирование таблицы
procedure Tab(n: Integer; var x, f:vec);
var
i: Integer;
begin
for i := 0 to n do
begin
Write('X', i:2, ',F',i:2,'?');
Readln(x[i], f[i]);
end;
end;

//Расчет коэффициентов сплайна
procedure CalcCoeff(n: Integer; x, f: vec; var c:vec);
var
i, j, m: Integer; a, b, r: Real; k: vec;
begin
k[1] := 0.0; c[1] := 0.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 CalcSpline(n: Integer; x, f, c:vec;
var p, p1, p2: Real; x1: Real);
var
i, j: Integer; a, b, d, q, r: Real;
begin
i := 1;
while (x1 > x[i]) and (i <> n) do
begin
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;
end;

begin
Write('n, x0, x9, h?');
Readln(n, x0, x9, h);
Tab(n, x, f);
CalcCoeff(n, x, f, c);
k := Round((x9 - x0) / h + 1.0);
x1 := x0;
for i := 0 to k do
begin
CalcSpline(n, x, f, c, p, p1, p2, x1);
Writeln(x1:5:5, ' ', p:5:5, ' ', p1:5:5, ' ', p2:5:5);
x1 := x1 + h;
end;
readln;
end.



5 comments:

  1. Вспоминаешь молодость?

    ReplyDelete
  2. Нет, Андрюх, я еще не бросил затею дописать диссер, а этот код мне понадобился для него :)

    ReplyDelete
  3. Доброго времени суток! Помогите пожалуйста с построением кубического сплайна. Задача передо мной стоит такая: Есть таблица значений X,Y и Z(x,y). допустим X=(0 1 2 3 4 5 6) Y=(4.1 2.4 3.0 4.3 3.6 5.2 5.9), а Z это пусть будет матрица любых значений 7*7. Программа должна находить в любой точке значение z (т.е. z(1.5,5)= ??). Это возможно? Пожалуйста наведи на верный путь..

    ReplyDelete
  4. Гляньте вот этот документ http://url4.ru/0cc

    ReplyDelete
  5. Большое спасибо за алгоритм. Сейчас сделаю в Excel.

    ReplyDelete