@@ -92,9 +92,9 @@ public NaturalSplineInterpolation(Vector<T> x, Vector<T> y, int degree = 3, IMat
9292 _degree = degree ;
9393 _decomposition = decomposition ;
9494 _numOps = MathHelper . GetNumericOperations < T > ( ) ;
95- _coefficients = new Vector < T > [ _degree ] ;
95+ _coefficients = new Vector < T > [ _degree + 1 ] ;
9696
97- for ( int i = 0 ; i < _degree ; i ++ )
97+ for ( int i = 0 ; i <= _degree ; i ++ )
9898 {
9999 _coefficients [ i ] = new Vector < T > ( x . Length - 1 ) ;
100100 }
@@ -125,7 +125,7 @@ public T Interpolate(T x)
125125 T dx = _numOps . Subtract ( x , _x [ i ] ) ;
126126 T result = _y [ i ] ;
127127
128- for ( int j = 1 ; j < _degree ; j ++ )
128+ for ( int j = 1 ; j <= _degree ; j ++ )
129129 {
130130 result = _numOps . Add ( result , _numOps . Multiply ( _coefficients [ j ] [ i ] , Power ( dx , j ) ) ) ;
131131 }
@@ -150,46 +150,69 @@ public T Interpolate(T x)
150150 private void CalculateCoefficients ( )
151151 {
152152 int n = _x . Length ;
153- Matrix < T > A = new Matrix < T > ( n , n ) ;
154- Vector < T > b = new Vector < T > ( n ) ;
155153
156- // Set up the system of equations
154+ // Calculate interval widths h[i] = x[i+1] - x[i]
155+ Vector < T > h = new Vector < T > ( n - 1 ) ;
157156 for ( int i = 0 ; i < n - 1 ; i ++ )
158157 {
159- T h = _numOps . Subtract ( _x [ i + 1 ] , _x [ i ] ) ;
160- A [ i , i ] = h ;
161- if ( i < n - 2 )
162- A [ i , i + 1 ] = _numOps . Multiply ( _numOps . FromDouble ( 2 ) , _numOps . Add ( h , _numOps . Subtract ( _x [ i + 2 ] , _x [ i + 1 ] ) ) ) ;
163- if ( i > 0 )
164- A [ i , i - 1 ] = h ;
165-
166- if ( i < n - 2 )
167- {
168- T dy1 = _numOps . Divide ( _numOps . Subtract ( _y [ i + 1 ] , _y [ i ] ) , h ) ;
169- T dy2 = _numOps . Divide ( _numOps . Subtract ( _y [ i + 2 ] , _y [ i + 1 ] ) , _numOps . Subtract ( _x [ i + 2 ] , _x [ i + 1 ] ) ) ;
170- b [ i ] = _numOps . Multiply ( _numOps . FromDouble ( 6 ) , _numOps . Subtract ( dy2 , dy1 ) ) ;
171- }
158+ h [ i ] = _numOps . Subtract ( _x [ i + 1 ] , _x [ i ] ) ;
172159 }
173160
174- // Apply natural spline conditions
161+ // Set up the tridiagonal system for natural cubic spline
162+ // We solve for second derivatives M[i] at each point
163+ Matrix < T > A = new Matrix < T > ( n , n ) ;
164+ Vector < T > b = new Vector < T > ( n ) ;
165+
166+ // Natural spline boundary conditions: M[0] = 0, M[n-1] = 0
167+ // Ensure boundary rows are completely set (all zeros except diagonal = 1)
168+ for ( int j = 0 ; j < n ; j ++ )
169+ {
170+ A [ 0 , j ] = _numOps . Zero ;
171+ A [ n - 1 , j ] = _numOps . Zero ;
172+ }
175173 A [ 0 , 0 ] = _numOps . One ;
176- A [ n - 1 , n - 1 ] = _numOps . One ;
177174 b [ 0 ] = _numOps . Zero ;
175+ A [ n - 1 , n - 1 ] = _numOps . One ;
178176 b [ n - 1 ] = _numOps . Zero ;
179177
180- // Solve the system
178+ // Interior equations (tridiagonal system)
179+ // h[i-1]*M[i-1] + 2*(h[i-1]+h[i])*M[i] + h[i]*M[i+1] = 6*((y[i+1]-y[i])/h[i] - (y[i]-y[i-1])/h[i-1])
180+ for ( int i = 1 ; i < n - 1 ; i ++ )
181+ {
182+ A [ i , i - 1 ] = h [ i - 1 ] ;
183+ A [ i , i ] = _numOps . Multiply ( _numOps . FromDouble ( 2 ) , _numOps . Add ( h [ i - 1 ] , h [ i ] ) ) ;
184+ A [ i , i + 1 ] = h [ i ] ;
185+
186+ T slope1 = _numOps . Divide ( _numOps . Subtract ( _y [ i ] , _y [ i - 1 ] ) , h [ i - 1 ] ) ;
187+ T slope2 = _numOps . Divide ( _numOps . Subtract ( _y [ i + 1 ] , _y [ i ] ) , h [ i ] ) ;
188+ b [ i ] = _numOps . Multiply ( _numOps . FromDouble ( 6 ) , _numOps . Subtract ( slope2 , slope1 ) ) ;
189+ }
190+
191+ // Solve the system for second derivatives M
181192 var decomposition = _decomposition ?? new LuDecomposition < T > ( A ) ;
182- Vector < T > m = MatrixSolutionHelper . SolveLinearSystem ( b , decomposition ) ;
193+ Vector < T > M = MatrixSolutionHelper . SolveLinearSystem ( b , decomposition ) ;
183194
184- // Calculate the coefficients
195+ // Calculate the spline coefficients for each segment
196+ // S(x) = a + b*(x-x[i]) + c*(x-x[i])^2 + d*(x-x[i])^3
185197 for ( int i = 0 ; i < n - 1 ; i ++ )
186198 {
187- T h = _numOps . Subtract ( _x [ i + 1 ] , _x [ i ] ) ;
199+ // a[i] = y [i]
188200 _coefficients [ 0 ] [ i ] = _y [ i ] ;
189- _coefficients [ 1 ] [ i ] = _numOps . Divide ( _numOps . Subtract ( _y [ i + 1 ] , _y [ i ] ) , h ) ;
190- _coefficients [ 1 ] [ i ] = _numOps . Subtract ( _coefficients [ 1 ] [ i ] , _numOps . Multiply ( _numOps . Divide ( h , _numOps . FromDouble ( 6 ) ) , _numOps . Add ( _numOps . Multiply ( _numOps . FromDouble ( 2 ) , m [ i ] ) , m [ i + 1 ] ) ) ) ;
191- _coefficients [ 2 ] [ i ] = _numOps . Divide ( m [ i ] , _numOps . FromDouble ( 2 ) ) ;
192- _coefficients [ 3 ] [ i ] = _numOps . Divide ( _numOps . Subtract ( m [ i + 1 ] , m [ i ] ) , _numOps . Multiply ( _numOps . FromDouble ( 6 ) , h ) ) ;
201+
202+ // b[i] = (y[i+1] - y[i])/h[i] - h[i]*(2*M[i] + M[i+1])/6
203+ T term1 = _numOps . Divide ( _numOps . Subtract ( _y [ i + 1 ] , _y [ i ] ) , h [ i ] ) ;
204+ T term2 = _numOps . Divide (
205+ _numOps . Multiply ( h [ i ] , _numOps . Add ( _numOps . Multiply ( _numOps . FromDouble ( 2 ) , M [ i ] ) , M [ i + 1 ] ) ) ,
206+ _numOps . FromDouble ( 6 ) ) ;
207+ _coefficients [ 1 ] [ i ] = _numOps . Subtract ( term1 , term2 ) ;
208+
209+ // c[i] = M[i]/2
210+ _coefficients [ 2 ] [ i ] = _numOps . Divide ( M [ i ] , _numOps . FromDouble ( 2 ) ) ;
211+
212+ // d[i] = (M[i+1] - M[i])/(6*h[i])
213+ _coefficients [ 3 ] [ i ] = _numOps . Divide (
214+ _numOps . Subtract ( M [ i + 1 ] , M [ i ] ) ,
215+ _numOps . Multiply ( _numOps . FromDouble ( 6 ) , h [ i ] ) ) ;
193216 }
194217 }
195218
0 commit comments