Need help with this matlab code for numerical methods. first time coding, I keep
ID: 3909468 • Letter: N
Question
Need help with this matlab code for numerical methods. first time coding, I keep getting error undefined function or Variable Xintx. If I change to Xint i get index exceeds array bounds. what should I do? how should I define the variable? Below is my code.
function Yint=CubicSplines(x,y,Xint)
n=length(x);
if(Xint(n))
error('out of x range!');
end
h=1:1:n;
b=1:1:n;
u=1:1:n;
v=1:1:n;
for i=1:n-1
h(i)=x(i+1)-x(i);
b(i)=y(i+1)-y(i)/h(i);
end
u(2)=2.0*(h(1)+h(2));
v(2)=6.0*(b(2)-b(1));
for i=3:n-1
u(i)=2.0*(h(i)+h(i-1))-h(i-1)^2/u(i-1);
v(i)=6.0*(b(i)-b(i-1))-h(i-1)*v(i-1)/u(i-1);
end
e(i)=0.0;
e(n)=0.0;
for i=n-1:-1:2
e(i)=(v(i)-h(i)*e(i+1))/u(i);
end
for i=1:n-1
if(Xint>=x(i)) && (Xint<=x(i+1))
m=i;
break
end
end
hm=x(m+1)-x(m);
xt1=Xint-x(m);
xt2=x(m+1)-Xint;
Yint=(e(m+1)*xt1^3+e(m)*xt2^3)/(6.0*hm)+xt1*(y(m+1)/hm-hm*e(m+1)/6.0)+xt2*(y(m)/hm-hm*e(m)/6.0);
end
Explanation / Answer
#include<stdio.h>
#include<math.h>
double F(double);
double CubicSpline(double *, double *, double *, double *, double *, double *, int, double);
void ThomasAlgorithm();
main()
{
double X[10],Y[10];
double A[10], B[10], C[10], D[10];
double A_Thomas[10], B_Thomas[10], C_Thomas[10], D_Thomas[10];
double Xint, Yint;
int N, i;
/* Set the number of points used in define the cubic spline. */
N = 4;
/* Assign the data points defining the cubic spline
interpolant. */
X[0] = 0.0;
Y[0] = F(X[0]);
X[1] = 0.5;
Y[1] = F(X[1]);
X[2] = 1.0;
Y[2] = F(X[2]);
X[3] = 1.5;
Y[3] = F(X[3]);
/* Specify the location to be interpolated. */
Xint = 0.85;
/* Determine the coefficients of the cubic spline
on each subinterval. This involves solving a
linear tridiagonal system of equations for the
coefficients of the second order term using the
Thomas algorithm. */
/* Step 1. Evaluate coefficient matrix of tridiagonal system. */
for (i=1; i<=N; i++) {
if (i == 1 || i == N) {
A_Thomas[i-1]=0.00;
D_Thomas[i-1]=1.00;
C_Thomas[i-1]=0.00;
B_Thomas[i-1]=0.00;
} else {
A_Thomas[i-1]=(X[i-1]-X[i-2])/3.00;
D_Thomas[i-1]=2.00*((X[i-1]-X[i-2])+(X[i]-X[i-1]))/3.00;
C_Thomas[i-1]=(X[i]-X[i-1])/3.0000;
B_Thomas[i-1]=(Y[i]-Y[i-1])/(X[i]-X[i-1])-
(Y[i-1]-Y[i-2])/(X[i-1]-X[i-2]);
}
}
/* Step 2. Solve tridiagonal system of equations using
subroutine THOMAS to determine coefficients C. */
ThomasAlgorithm(A_Thomas, D_Thomas, C_Thomas, B_Thomas, C, N);
/* Step 3. Determine the other coefficients A, B, and D. */
for (i=1; i<=N-1; i++) {
A[i-1]=Y[i-1];
B[i-1]=(Y[i]-Y[i-1])/(X[i]-X[i-1])-
(2.00*C[i-1]+C[i])*(X[i]-X[i-1])/3.00;
D[i-1]=(C[i]-C[i-1])/(3.00*(X[i]-X[i-1]));
}
/* Perform the cubic spline interpolation. */
Yint = CubicSpline(X, Y, A, B, C, D, N, Xint);
/* Output the results. */
printf("Results of cubic spline interpolation: ");
printf("x = %11.8f y = %11.8f ", Xint, Yint);
/* End of program. */
return 0;
}
/* Subroutine that defines the function to be interpolated. */
double F(double x)
{
double f;
f = sin(x);
return f;
}
/* Subroutine that performs the cubic spline interpolation. */
double CubicSpline(double *X, double *Y, double *A , double *B, double *C, double *D, int N, double Xint)
{
double y;
int interval;
/* Check to see if XT is within the valid range
for the cubic spline. */
if (Xint < X[0] || Xint > X[N-1]) {
y = 0.00;
return y;
}
/* Determine which subinterval contains Xint. */
interval = 1;
while (!(Xint >= X[interval-1] && Xint < X[interval]) && interval < N-1) {
++interval;
}
/* Evaluate the cubic polynomial fit to the data for
subinterval containing the location of interest. */
y = A[interval-1] +
B[interval-1]*(Xint-X[interval-1]) +
C[interval-1]*pow(Xint-X[interval-1], 2) +
D[interval-1]*pow(Xint-X[interval-1], 3);
/* Return the interpolated value. */
return y;
}
/* Subroutine applies the Thomas algorithm to solve a tridiagonal
linear system of equations. */
void ThomasAlgorithm(double *A, double *D, double *C, double *B, double *X, int N)
{
int i;
/* Perform forward elimination. */
D[0]=D[0];
B[0]=B[0];
for (i=2; i<=N; i++) {
D[i-1]=D[i-1]-A[i-1]*C[i-2]/D[i-2];
B[i-1]=B[i-1]-A[i-1]*B[i-2]/D[i-2];
}
/* Perform back substitution. */
X[N-1]=B[N-1]/D[N-1];
for (i=N-1; i>=1; i--) {
X[i-1]=(B[i-1]-C[i-1]*X[i])/D[i-1];
}
}
This code will solve index out of bound problem.
Thank You :)