具體推導看書<<數值分析>>
code：
#include <iostream>
using namespace std;
const int MAXN = 100;
int n;
double x[MAXN], y[MAXN]; //下標從0..n
double alph[MAXN], beta[MAXN], a[MAXN], b[MAXN];
double h[MAXN];
double m[MAXN]; //各點的一階導數;
inline double sqr(double pa) 
{
return pa * pa;
}
double sunc(double p, int i) {
return (1 + 2 * (p - x[i]) / (x[i + 1] - x[i])) * sqr((p - x[i + 1]) / (x[i + 1] - x[i])) * y[i]
+ (1 + 2 * (p - x[i + 1]) / (x[i] - x[i + 1])) * sqr((p - x[i]) / (x[i + 1] - x[i])) * y[i + 1]
+ (p - x[i]) * sqr((p - x[i + 1]) / (x[i] - x[i + 1])) * m[i]
+ (p - x[i + 1]) * sqr((p - x[i]) / (x[i + 1] - x[i])) * m[i + 1];
}
int main() {
int i, j;
double xx;
freopen("threeInsert.in", "r", stdin);
scanf("%d", &n);
for (i = 0; i <= n; i++) scanf("%lf%lf", &x[i], &y[i]);
// scanf("%lf%lf", &m[0], &m[n]);
for (i = 0; i <= n - 1; i++) h[i] = x[i + 1] - x[i];
//第一種邊界條件
//alph[0] = 0; alph[n] = 1; beta[0] = 2 * m[0]; beta[n] = 2 * m[n];
//第二種邊界條件
alph[0] = 1; alph[n] = 0; beta[0] = 3 * (y[1] - y[0]) / h[0]; beta[n] = 3 * (y[n] - y[n - 1] / h[n - 1]);
for (i = 1; i <= n - 1; i++) {
alph[i] = h[i - 1] / (h[i - 1] + h[i]);
beta[i] = 3 * ((1 - alph[i]) * (y[i] - y[i - 1]) / h[i - 1] + alph[i] * (y[i + 1] - y[i]) / h[i]);
}
a[0] = - alph[0] / 2; b[0] = beta[0] / 2;
for (i = 1; i <= n; i++) {
a[i] = - alph[i] / (2 + (1 - alph[i]) * a[i - 1]);
b[i] = (beta[i] - (1 - alph[i]) * b[i - 1]) / (2 + (1 - alph[i]) * a[i - 1]);
}
m[n + 1] = 0;
for (i = n; i >= 0; i--) {
m[i] = a[i] * m[i + 1] + b[i];
}
scanf("%lf", &xx);
for (i = 0; i < n; i++) {
if (xx >= x[i] && xx <= x[i + 1]) break;
}
printf("%lf\n", sunc(xx, i));
return 0;
} posted on 2007-10-20 13:07
豪 閱讀(3527)
評論(4) 編輯 收藏 引用 所屬分類:
算法&ACM