Полином Ньютона на Паскале и C++
Предупреждение: полиномы Ньютона (их два) - вещь гнилая. 1-й работает плохо в конце отрезка, 2-й - в начале. На практике вряд ли используются
Письмо:
Добрый день! Я являюсь постоянным читателем http://nickolay.info. У меня к вам огромная просьба, не могли бы вы помочь мне разобраться с задачей на Pascal связанная с численными методами, читал у вас "Реализация основных методов интерполяции функции одной переменной", но это не совсем то что мне нужно. Если бы вы помогли разобраться со следующей задачей я был бы вам очень признателен: Построить интерполяционный многочлен Ньютона по известным в узлах х0,х1,...,х5 значениями f(x). Вычислить значения этого многочлена в равномерной сетке точек [a,b] с шагом h и сравнить их со значением функции f(x) в этих точках. Использовать схему интерполяции вперед.
Функция f(x)=x+exp(-x^2).
Значения функции:
Х0=0,2
Х1=0,5
Х2=1
Х3=1,3
Х4=1,5
Х5=1.9
Отрезок:
a=0
b=2
Шаг: h=0,05
Ответ: (вдруг кому пригодится, всё ж минуты 4 "труда" :)
Добрый день, насколько я ориентируюсь на своём сайте, полином Ньютона на Паскале есть вот тут: http://nickolay.info/study/pas_methods.html оттуда и можно надёргать кода, как-то вот так, не проверял толком:
Uses Crt; Type vector=array[1..50] of real; Procedure Interpol (var X,Y:vector;n,nom:integer;a:real;var f,r:real); Var D1,D2,D3,D4:vector; h,q:real; k,i:integer; Begin for i:=1 to n-1 do d1[i]:=y[i+1]-y[i]; for i:=1 to n-2 do d2[i]:=d1[i+1]-d1[i]; for i:=1 to n-3 do d3[i]:=d2[i+1]-d2[i]; for i:=1 to n-4 do d4[i]:=d3[i+1]-d3[i]; h:=abs(x[2]-x[1]); case nom of 1:begin i:=1; while a>x[i] do i:=i+1; k:=i-1; q:=(a-x[k])/h; f:=y[k]+d1[k]*q+d2[k]*q*(q-1)/2+ d3[k]*q*(q-1)*(q-2)/6; R:=abs(d4[k]*q*(q-1)*(q-2)*(q-3)/24); end; 2:begin i:=1; while a>x[i] do i:=i+1; k:=i; q:=(a-x[k])/h; f:=y[k]+d1[k-1]*q+d2[k-2]*q*(q+1)/2+ d3[k-3]*q*(q+1)*(q+2)/6; R:=abs(d4[k-4]*q*(q+1)*(q+2)*(q+3)/24); end; end; End; function f(x:real):real; begin f:=x+exp(-sqr(x)); end; Var X,Y:vector; n,nom,i:integer; a,b,h,xx,ff,R:real; Begin ClrScr; write('Введите размер таблицы n='); readln(n); writeln('Введите массив X='); for i:=1 to n do read(x[i]); write('массив Y='); for i:=1 to n do begin y[i]:=f(x[i]); write (y[i]:18:6); end; writeln; write('Введите интервал [a,b]='); readln(a,b); write('Введите шаг h='); readln(h); write('Введите номер интерполяционной ', 'формулы Ньютона: nom='); readln(nom); xx:=a; writeln ('X':18,'Вычислено':18,'Точное':18,'остаток':18); while xx<=b do begin Interpol(X,Y,n,nom,xx,ff,R); writeln(xx:18:6,ff:18:6,f(xx):18:6,R:18:6); xx:=xx+h; end; reset (input); readln; End.
Скриншот тестового прогона на этих данных (Турбо Паскаль 7.1 в консольке):
Скриншот Полином Ньютона на Паскале
Узлы интерполяции по оси X предполагаются равноотстоящими. Ясно, что вне интервала ]Xmin,Xmax[
вычисленные значения смысла не имеют. Также может и глюк где.
Вот C++-реализация первой интерполяционный формулы Ньютона (для левой половины интервала), запускалась в консоли Visual Studio.
#include <iostream> using namespace std; struct Data { double x, y; }; double u_cal (double u, int n) { double temp = u; for (int i = 1; i < n; i++) temp *= (u - i); return temp; } double fact(int n) { double f = 1.; for (int i = 2; i <= n; i++) f *= i; return f; } int main() { //Количество исходных точек: const int n = 4; //Исходные значения (x,y): Data f[n] = { { 0, 2 }, { 1, 3 }, { 2, 12 }, { 5, 147 } }; //Таблица конечных разностей (вперёд) double y[n][n]; for (int i = 0; i < n; i++) y[i][0]=f[i].y; for (int i = 1; i < n; i++) { for (int j = 0; j < n - i; j++) y[j][i] = y[j + 1][i - 1] - y[j][i - 1]; } //Вывод таблицы значений от x1 до x2 с шагом dx: double x1 = 0, x2 = 5, dx = 0.125; cout.width(20); cout << "X"; cout.width(20); cout << "Y" << endl; for (double x = x1; x <= x2; x += dx) { double sum = y[0][0]; double u = (x - f[0].x) / (f[1].x - f[0].x); for (int i = 1; i < n; i++) { sum += (u_cal(u, i) * y[0][i]) / fact(i); } cout.width(20); cout.precision(14); cout << x; cout.width(20); cout.precision(14); cout << sum << endl; } cin.get(); return 0; }
Вторая интерполяционная формула Ньютона (для правой половины интервала).
#include <iostream> using namespace std; struct Data { double x, y; }; double u_cal (double u, int n) { double temp = u; for (int i = 1; i < n; i++) temp *= (u + i); return temp; } double fact(int n) { double f = 1.; for (int i = 2; i <= n; i++) f *= i; return f; } int main() { //Количество исходных точек: const int n = 4; //Исходные значения (x,y): Data f[n] = { { 0, 2 }, { 1, 3 }, { 2, 12 }, { 5, 147 } }; //Таблица конечных разностей (назад) double y[n][n]; for (int i = 0; i < n; i++) y[i][0]=f[i].y; for (int i = 1; i < n; i++) { for (int j = n - 1; j >= i; j--) y[j][i] = y[j][i - 1] - y[j - 1][i - 1]; } //Вывод таблицы значений от x1 до x2 с шагом dx: double x1 = 0, x2 = 5, dx = 0.125; cout.width(20); cout << "X"; cout.width(20); cout << "Y" << endl; for (double x = x1; x <= x2; x += dx) { double sum = y[n-1][0]; double u = (x - f[n-1].x) / (f[1].x - f[0].x); for (int i = 1; i < n; i++) { sum += (u_cal(u, i) * y[n-1][i]) / fact(i); } cout.width(20); cout.precision(14); cout << x; cout.width(20); cout.precision(14); cout << sum << endl; } cin.get(); return 0; }
19.03.2013, 16:53 [32783 просмотра]