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

@MayaNash · Регистрация: 24.08.2011

Студворк — интернет-сервис помощи студентам

Добрый вечер.
В общем, тему не знаю, а программу написать надо. Нашла статейку, вроде понятно (см. вложение, взято из http://bookzooka.com/book/753-... ajnov.html ). Но дан пример процедуры:
Входные параметры: a,b – интервал дифференцирования; x – точка дифференцирования; n – число точек дифференцирования; fun – вид функции.
Выходные параметры: y_1,y_2 – значения первой и второй производной в точке x соответственно.
Продифференцировать функцию $https://www.cyberforum.ru/cgi-bin/latex.cgi?{e}^{{x}^{2}}$ , количество точек разбиения 5, значение точки х=1.

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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;

Формулы в процедуре неправильные, или я какой-то не тот способ почитала? Или это формулы интегрирования, что вроде бы следует из названия? Дано точно в примере процедуры дифференцирования.

@S_el · 13.01.2015, 20:45

Так а в чем вопрос?Соответствует ли программа,описанным формулам?

@MayaNash · 13.01.2015, 20:53 **[ТС]**

Имеют ли формулы в программе вообще отношение к сплайновому дифференцированию? А то меня смущает такая разница. Может я не тот способ нашла.

@S_el · 13.01.2015, 21:12

Сообщение от Керра

Или это формулы интегрирования, что вроде бы следует из названия?

Формулы для дифференцирования.

Сообщение от Керра

Имеют ли формулы в программе вообще отношение к сплайновому дифференцированию? А то меня смущает такая разница. Может я не тот способ нашла.

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

@MayaNash · 14.01.2015, 12:05 **[ТС]**

Попробовала таки те формулы. Первая производная считается верно, проверила на функциях arctg(1/x) и (x*x+2)/3. Но вторая считается неправильно, даже если ее считать через формулу первой производной, запомнив значения первой производной, т.е. находить производную от первой производной. А по той формуле второй производной получаются большие числа из-за деления последнего слагаемого на слишком маленькое число.

Кликните здесь для просмотра всего текста

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
uses crt;
 
type
  p_function = function(x: real): real;
 
//----------------------------------------------------------
function f1(x: real): real;
begin
  f1 := arctan(1/x);
end;
//----------------------------------------------------------
procedure splines(f: p_function; a, b, x0: real; n: integer;
                  var dy1_, dy2_: real);
var
  x, y, dy, dy1: array of real;
  h: real;
  i: integer;
begin
  h := (b - a)/n;
  setlength(x, n-1);
  setlength(y, n-1);
  setlength(dy, n-1);
  setlength(dy1, n-1);
  x[0] := a;
  y[0] := f(a);
  dy[0] := 0;
  for i := 1 to n-2 do
  begin
    x[i] := x[i-1]+h;
    y[i] := f(x[i]);
    dy[i] := (y[i] - y[i-1])/h;
    dy1[i] := (dy[i] - dy[i-1])/h;
  end;
 
  i := trunc((x0 - a)/h + 0.5);
  dy1_ := dy[i-1]*(x[i] -x0)*(2*x[i-1] + x[i] - 3*x0)/sqr(h)
       - dy[i]*(x0 - x[i-1])*(2*x[i] + x[i-1] - 3*x0)/sqr(h)
       + 6*(y[i] - y[i-1])*(x[i] - x0)*(x0 - x[i-1])/(h*h*h);
  dy2_ := dy1[i-1]*(x[i] -x0)*(2*x[i-1] + x[i] - 3*x0)/sqr(h)
       - dy1[i]*(x0 - x[i-1])*(2*x[i] + x[i-1] - 3*x0)/sqr(h)
       + 6*(y[i] - y[i-1])*(x[i] - x0)*(x0 - x[i-1])/(h*h*h);
  {dy2_ := -2*dy[i-1]*(2*x[i] + x[i-1] - 3*x0)/sqr(h)
       - 2*dy[i]*(2*x[i-1] + x[i] - 3*x0)/sqr(h)
       + 6*(y[i] - y[i-1])*(x[i] - x[i-1] - 2*x0)/sqr(h);
  writeln(i);
  writeln(y[i]:0:10);
  writeln(y[i-1]:0:10);
  writeln(x[i]:0:3);
  writeln(x[i-1]:0:3);
  writeln(h:0:3);
  writeln(-2*dy[i-1]*(2*x[i] + x[i-1] - 3*x0)/sqr(h):0:3);
  writeln(2*dy[i]*(2*x[i-1] + x[i] - 3*x0)/sqr(h):0:3);
  writeln(6*(y[i] - y[i-1])*(x[i] - x[i-1] - 2*x0)/sqr(h):0:3);
  writeln('---------------');}
end;
//----------------------------------------------------------
 
var
  dy1, dy2: real;
 
begin
  clrscr;
  splines(@f1, 1, 4, 3, 400, dy1, dy2);
  writeln(dy1:0:5);
  writeln(dy2:0:5);
  readkey;
end.

Почему не получается?

Ilot · 14.01.2015, 12:48

Ошибка здесь:
dy2_ := -2*dy[i-1]*(2*x[i] + x[i-1] - 3*x0)/sqr(h)
- 2*dy[i]*(2*x[i-1] + x[i] - 3*x0)/sqr(h)
+ 6*(y[i] - y[i-1])*(x[i] + x[i-1] - 2*x0)/sqr(h);
Это ошибка на сайте.

@MayaNash · 14.01.2015, 12:59 **[ТС]**

Ilot, $https://www.cyberforum.ru/cgi-bin/latex.cgi?{f}_{1}=arctg\frac{1}{x}$
$https://www.cyberforum.ru/cgi-bin/latex.cgi?{f}_{1}'=-\frac{1}{{x}^{2}+1}$
$https://www.cyberforum.ru/cgi-bin/latex.cgi?{f}_{1}''=\frac{2x}{{({x}^{2}+1)}^{2}}$
$https://www.cyberforum.ru/cgi-bin/latex.cgi?{f}_{1}'(3)=-0.1$ - программа отвечает так же
Но $https://www.cyberforum.ru/cgi-bin/latex.cgi?{f}_{1}''(3)=0.6$ , а программа отвечает -26.48656 по вашей формуле. По той было гораздо хуже, но все равно не то выходит

Ilot · 14.01.2015, 13:05

Керра, ну во-первых, производная арктангенса без минуса. Во-вторых, там, на сайте, еще одна ошибка. В знаменателе в последней дроби должна стоять третья степень h иначе опять получим не соответствие размерностей.

@MayaNash · 14.01.2015, 13:12 **[ТС]**

Ilot, но там 1/x: $https://www.cyberforum.ru/cgi-bin/latex.cgi?(arctg\frac{1}{x})'=\frac{1}{1+\frac{1}{{x}^{2}}}*(-\frac{1}{{x}^{2}})=-\frac{1}{{x}_{2}+1}$

Ilot · 14.01.2015, 13:14

Керра, виноват моя госпажа.

Не досмотрел.

@MayaNash · 14.01.2015, 13:18 **[ТС]**

Ilot,

Вопрос в том что теперь делать... Третий день это пытаюсь отладить...

Ilot · 14.01.2015, 13:20

Керра, ну а степень то исправили?
Хм... сейчас сам попробую запилить код. Посмотрим, что там...

@MayaNash · 14.01.2015, 13:26 **[ТС]**

Теперь -0

Кликните здесь для просмотра всего текста

Pascal
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
uses crt;
 
type
  p_function = function(x: real): real;
 
//----------------------------------------------------------
function f1(x: real): real;
begin
  f1 := arctan(1/x);
end;
//----------------------------------------------------------
procedure splines(f: p_function; a, b, x0: real; n: integer;
                  var dy1_, dy2_: real);
var
  x, y, dy, dy1: array of real;
  h: real;
  i: integer;
begin
  h := (b - a)/n;
  setlength(x, n-1);
  setlength(y, n-1);
  setlength(dy, n-1);
  setlength(dy1, n-1);
  x[0] := a;
  y[0] := f(a);
  dy[0] := 0;
  for i := 1 to n-2 do
  begin
    x[i] := x[i-1]+h;
    y[i] := f(x[i]);
    dy[i] := (y[i] - y[i-1])/h;
    //dy1[i] := (dy[i] - dy[i-1])/h;
  end;
 
  i := trunc((x0 - a)/h + 0.5);
  dy1_ := dy[i-1]*(x[i] -x0)*(2*x[i-1] + x[i] - 3*x0)/sqr(h)
       - dy[i]*(x0 - x[i-1])*(2*x[i] + x[i-1] - 3*x0)/sqr(h)
       + 6*(y[i] - y[i-1])*(x[i] - x0)*(x0 - x[i-1])/(h*h*h);
  {dy2_ := dy1[i-1]*(x[i] -x0)*(2*x[i-1] + x[i] - 3*x0)/sqr(h)
       - dy1[i]*(x0 - x[i-1])*(2*x[i] + x[i-1] - 3*x0)/sqr(h)
       + 6*(y[i] - y[i-1])*(x[i] - x0)*(x0 - x[i-1])/(h*h*h);}
  dy2_ := -2*dy[i-1]*(2*x[i] + x[i-1] - 3*x0)/sqr(h)
       - 2*dy[i]*(2*x[i-1] + x[i] - 3*x0)/sqr(h)
       + 6*(y[i] - y[i-1])*(x[i] + x[i-1] - 2*x0)/(h*h*h);
  {writeln(i);
  writeln(y[i]:0:10);
  writeln(y[i-1]:0:10);
  writeln(x[i]:0:3);
  writeln(x[i-1]:0:3);
  writeln(h:0:3);
  writeln(-2*dy[i-1]*(2*x[i] + x[i-1] - 3*x0)/sqr(h):0:3);
  writeln(2*dy[i]*(2*x[i-1] + x[i] - 3*x0)/sqr(h):0:3);
  writeln(6*(y[i] - y[i-1])*(x[i] - x[i-1] - 2*x0)/sqr(h):0:3);
  writeln('---------------');}
end;
//----------------------------------------------------------
 
var
  dy1, dy2: real;
 
begin
  clrscr;
  splines(@f1, 1, 4, 3, 400, dy1, dy2);
  writeln(dy1:0:5);
  writeln(dy2:0:10);
  readkey;
end.

Ilot · 14.01.2015, 13:52

Керра, все там в порядке. Вот код на с++. Меняйте к от 0 до n. Совпадение вплоть до точности представления:

Кликните здесь для просмотра всего текста

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#include <iostream>
#include <cmath>
int main() {
    const int n = 20;
    double y[n + 1];
    double m[n + 1];
    double a = 0.0;
    double b = 1.0;
    double h = (b - a) / n;
    for(int i = 0; i < n + 1; ++i) {
        y[i] = sin(a + i * h);
        m[i] = cos(a + i * h);
    }
    int k = 3;
    double x = a + 0.7 * h + (k - 1) * h;
    double dy = (
                    - 2 * m[k-1] * (3 * a + (3*k-1) * h - 3*x)
                    - 2 * m[k]   * (3 * a + (3*k-2) * h - 3*x)
                ) /(h*h) + 6 * (y[k] - y[k-1]) * (2*a + (2*k-1) * h - 2 * x) / (h * h * h);
    std::cout << "sin''(x) : " << -sin(x) << std::endl;
    std::cout << "phi(x)   : " << dy << std::endl;
 
    return 0;
}

Должно быть у вас где-то закралась ошибка.

@MayaNash · 14.01.2015, 14:37 **[ТС]**

Ilot, а если без формулы производной? Вы вычисляете m[i] с помощью формулы а по идее эта формула неизвестна, и функция может быть абсолютно любой, в том числе заданной значениями функции в точках.

Ilot · 14.01.2015, 15:26

Керра, ошибка округления. Дело в том, что в данной формуле получаются разности очень близких чисел. Отсюда и такой результат. Вообще-то далее в статье указанно как находить коэффициенты m. Явно простая формула через конечные разности здесь не применима.

@MayaNash · 14.01.2015, 15:57 **[ТС]**

А я надеялась их избежать...

@vioup · 22.02.2017, 00:35

Вы в конечном итоге решили эту задачу?

@MayaNash · 22.02.2017, 00:43 **[ТС]**

vioup, да как бы не помню уже

Боже, какой у меня был ужасный стиль написания кода... Смотреть больно...

@vioup · 22.02.2017, 09:45

Понятно. Сейчас сам страдаю от этой "процедуры". Уже и писал преподавателю в СДО, полное игнорирование. Сколько лет прошло, а нормальную процедуру написать было для них не вариант.

Новые блоги и статьи Все статьи Все блоги /
Оттенки серого Argus19 18.03.2026 Оттенки серого Нашёл в интернете 3 прекрасных модуля: Модуль класса открытия диалога открытия/ сохранения файла на Win32 API; Модуль класса быстрого перекодирования цветного изображения в оттенки. . .	SDL3 для Desktop (MinGW): Рисуем цветные прямоугольники с помощью рисовальщика SDL3 на Си и C++ 8Observer8 17.03.2026 Содержание блога Финальные проекты на Си и на C++: finish-rectangles-sdl3-c. zip finish-rectangles-sdl3-cpp. zip	Символические и жёсткие ссылки в Linux. algri14 15.03.2026 Существует два типа ссылок — символические и жёсткие. Ссылка в Linux — это запись в каталоге, которая может указывать либо на inode «файла-ИСТОЧНИКА», тогда это будет «жёсткая ссылка» (hard link),. . .	[Owen Logic] Поддержание уровня воды в резервуаре количеством включённых насосов: моделирование и выбор регулятора ФедосеевПавел 14.03.2026 Поддержание уровня воды в резервуаре количеством включённых насосов: моделирование и выбор регулятора ВВЕДЕНИЕ Выполняя задание на управление насосной группой заполнения резервуара,. . .
делаю науч статью по влиянию грибов на сукцессию anaschu 13.03.2026 прикрепляю статью	SDL3 для Desktop (MinGW): Создаём пустое окно с нуля для 2D-графики на SDL3, Си и C++ 8Observer8 10.03.2026 Содержание блога Финальные проекты на Си и на C++: hello-sdl3-c. zip hello-sdl3-cpp. zip Результат:	Установка CMake и MinGW 13.1 для сборки С и C++ приложений из консоли и из Qt Creator в EXE 8Observer8 10.03.2026 Содержание блога MinGW - это коллекция инструментов для сборки приложений в EXE. CMake - это система сборки приложений. Здесь описаны базовые шаги для старта программирования с помощью CMake и. . .	Как дизайн сайта влияет на конверсию: 7 решений, которые реально повышают заявки Neotwalker 08.03.2026 Многие до сих пор воспринимают дизайн сайта как “красивую оболочку”. На практике всё иначе: дизайн напрямую влияет на то, оставит человек заявку или уйдёт через несколько секунд. Даже если у вас. . .

@S_el 2444 / 1842 / 406 Регистрация: 15.12.2013 Сообщений: 8,243
	13.01.2015, 20:45
	Так а в чем вопрос?Соответствует ли программа,описанным формулам? 0

@MayaNash 1296 / 470 / 151 Регистрация: 24.08.2011 Сообщений: 2,249
	13.01.2015, 20:53 [ТС]
	Имеют ли формулы в программе вообще отношение к сплайновому дифференцированию? А то меня смущает такая разница. Может я не тот способ нашла. 0

Ilot $Эксперт по математике/физике$ 2224 / 1426 / 420 Регистрация: 16.05.2013 Сообщений: 3,647 Записей в блоге: 6
	14.01.2015, 12:48
	Ошибка здесь: dy2_ := -2dy[i-1](2x[i] + x[i-1] - 3x0)/sqr(h) - 2dy[i](2x[i-1] + x[i] - 3x0)/sqr(h) + 6(y[i] - y[i-1])(x[i] + x[i-1] - 2*x0)/sqr(h); Это ошибка на сайте. 0

@MayaNash 1296 / 470 / 151 Регистрация: 24.08.2011 Сообщений: 2,249
	14.01.2015, 12:59 [ТС]
	Ilot, $https://www.cyberforum.ru/cgi-bin/latex.cgi?{f}_{1}=arctg\frac{1}{x}$ $https://www.cyberforum.ru/cgi-bin/latex.cgi?{f}_{1}'=-\frac{1}{{x}^{2}+1}$ $https://www.cyberforum.ru/cgi-bin/latex.cgi?{f}_{1}''=\frac{2x}{{({x}^{2}+1)}^{2}}$ $https://www.cyberforum.ru/cgi-bin/latex.cgi?{f}_{1}'(3)=-0.1$ - программа отвечает так же Но $https://www.cyberforum.ru/cgi-bin/latex.cgi?{f}_{1}''(3)=0.6$ , а программа отвечает -26.48656 по вашей формуле. По той было гораздо хуже, но все равно не то выходит 0

Ilot $Эксперт по математике/физике$ 2224 / 1426 / 420 Регистрация: 16.05.2013 Сообщений: 3,647 Записей в блоге: 6
	14.01.2015, 13:05
	Керра, ну во-первых, производная арктангенса без минуса. Во-вторых, там, на сайте, еще одна ошибка. В знаменателе в последней дроби должна стоять третья степень h иначе опять получим не соответствие размерностей. 0

@MayaNash 1296 / 470 / 151 Регистрация: 24.08.2011 Сообщений: 2,249
	14.01.2015, 13:12 [ТС]
	Ilot, но там 1/x: $https://www.cyberforum.ru/cgi-bin/latex.cgi?(arctg\frac{1}{x})'=\frac{1}{1+\frac{1}{{x}^{2}}}*(-\frac{1}{{x}^{2}})=-\frac{1}{{x}_{2}+1}$ 0

Ilot $Эксперт по математике/физике$ 2224 / 1426 / 420 Регистрация: 16.05.2013 Сообщений: 3,647 Записей в блоге: 6
	14.01.2015, 13:14
	Керра, виноват моя госпажа. Не досмотрел. 0

@MayaNash 1296 / 470 / 151 Регистрация: 24.08.2011 Сообщений: 2,249
	14.01.2015, 13:18 [ТС]
	Ilot, Вопрос в том что теперь делать... Третий день это пытаюсь отладить... 0

Ilot $Эксперт по математике/физике$ 2224 / 1426 / 420 Регистрация: 16.05.2013 Сообщений: 3,647 Записей в блоге: 6
	14.01.2015, 13:20
	Керра, ну а степень то исправили? Хм... сейчас сам попробую запилить код. Посмотрим, что там... 0

@MayaNash 1296 / 470 / 151 Регистрация: 24.08.2011 Сообщений: 2,249
	14.01.2015, 14:37 [ТС]
	Ilot, а если без формулы производной? Вы вычисляете m[i] с помощью формулы а по идее эта формула неизвестна, и функция может быть абсолютно любой, в том числе заданной значениями функции в точках. 0

Ilot $Эксперт по математике/физике$ 2224 / 1426 / 420 Регистрация: 16.05.2013 Сообщений: 3,647 Записей в блоге: 6
	14.01.2015, 15:26
	Керра, ошибка округления. Дело в том, что в данной формуле получаются разности очень близких чисел. Отсюда и такой результат. Вообще-то далее в статье указанно как находить коэффициенты m. Явно простая формула через конечные разности здесь не применима. 0

@MayaNash 1296 / 470 / 151 Регистрация: 24.08.2011 Сообщений: 2,249
	14.01.2015, 15:57 [ТС]
	А я надеялась их избежать... 0

@vioup 0 / 0 / 0 Регистрация: 07.04.2015 Сообщений: 55
	22.02.2017, 00:35
	Вы в конечном итоге решили эту задачу? 0

@MayaNash 1296 / 470 / 151 Регистрация: 24.08.2011 Сообщений: 2,249
	22.02.2017, 00:43 [ТС]
	vioup, да как бы не помню уже Боже, какой у меня был ужасный стиль написания кода... Смотреть больно... 0

@vioup 0 / 0 / 0 Регистрация: 07.04.2015 Сообщений: 55
	22.02.2017, 09:45
	Понятно. Сейчас сам страдаю от этой "процедуры". Уже и писал преподавателю в СДО, полное игнорирование. Сколько лет прошло, а нормальную процедуру написать было для них не вариант. 0