1. // interpolacja_newtona_11k2.cpp : Defines the entry point for the console application.
2. //
3.  
4. #include "stdafx.h"
5. #include <cmath>
6. #include <iostream>
7.  
8. using namespace std;
9.  
10. double f(double x);
11. double W(double x1, double x2, double x3, double x, double a0, double a1, double a2, double a3);
12.  
13. int _tmain(int argc, _TCHAR* argv[])
14. {
15.         double x1, x2, x3, dx, er, ermax, erdop, xStart, xEnd, xp, a0, a1, a2, a3;
16.         int m = 3; // przedzialy
17.         double d[5][5];
18.         double x[5];
19.  
20.         xStart = 0.0;
21.         xEnd = 2.0*asin(1.0);
22.  
23.         cout << endl << "Podaj blad=";
24.         cin >> erdop;
25.  
26.         do
27.         {
28.                 ermax = 0.0;
29.                 xStart = 0.0;
30.                 xEnd = 2.0*asin(1.0);
31.                 dx = (xEnd - xStart) / m;
32.  
33.                 x[1] = xStart;
34.                 x[4] = xEnd / m;
35.                 x[2] = x[1] + (x[4] - x[1]) / 3;
36.                 x[3] = x[1] + (x[4] - x[1]) * 2 / 3;
37.                 xp = x[1] + dx / 2.0;
38.                
39.  
40.                 for (int i = 0; i<m; i++)
41.                 {
42.                         d[1][1] = f(x[1]);
43.                         d[2][1] = f(x[2]);
44.                         d[3][1] = f(x[3]);
45.                         d[4][1] = f(x[4]);
46.                        
47.                         for (int i = 2; i < 5; i++)
48.                         {
49.                                 for (int j = 2; j < 5; j++)
50.                                 {
51.                                         d[i][j] = (d[i][j - 1] - d[i - 1][j - 1]) / (x[i] - x[i - j + 1]);
52.                                 }
53.                         }
54.  
55.                         a0 = d[1][1];
56.                         a1 = d[2][2];
57.                         a2 = d[3][3];
58.                         a3 = d[4][4];
59.  
60.                         er = f(xp) - W(x[1], x[2], x[3], xp, a0, a1, a2, a3);
61.  
62.                         if (er<0)
63.                                 er = -er;
64.  
65.                         if (er>ermax)
66.                                 ermax = er;
67.  
68.                         x[1] = x[4];
69.                         x[4] = x[1] + dx;
70.                         x[2] = x[1] + (x[4] - x[1]) / 3;
71.                         x[3] = x[1] + (x[4] - x[1]) * 2 / 3;
72.                         xp = x[1] + dx / 2.0;
73.                 }
74.                 m++;
75.         } while (ermax>erdop);
76.  
77.         printf("\nliczba przedzialow=%d blad=%f", m - 1, ermax);
78.         cout << endl;
79.         system("pause");
80.         return 0;
81. }
82.  
83. double f(double x)
84. {
85.         return sin(x);
86. }
87. //---------------------------------------
88. double W(double x1, double x2, double x3, double x, double a0, double a1, double a2, double a3)
89. {
90.         return  a0+ a1*(x - x1) + a2*(x - x1)*(x - x2) + a3*(x - x1)*(x - x2)*(x - x3);
91. }