三次样条插值的C程序(很全啊)

三次样条插值 C C++ 程序

三次样条插值C/C++程序（自己整理的）

具体推导看书<<数值分析>>

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;

三次样条插值 C C++ 程序

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));

三次样条插值 C C++ 程序

return 0;

}

#include<iostream>

#include<iomanip>

using namespace std;

const int MAX = 50;

float x[MAX], y[MAX], h[MAX];//变量设置：x为各点横坐标；y为各点纵坐标；h为步长

float c[MAX], a[MAX], fxym[MAX];

float f(int x1, int x2, int x3)/*****************求差分函数(含三个参数)****************************/

{

float a = (y[x3] - y[x2]) / (x[x3] - x[x2]);

float b = (y[x2] - y[x1]) / (x[x2] - x[x1]);

return (a - b)/(x[x3] - x[x1]);

}

void cal_m(int n)/***********************用追赶法求解出弯矩向量M……***************************/

{

float B[MAX];

B[0] = c[0] / 2;

for(int i = 1; i < n; i++)

B[i] = c[i] / (2 - a[i]*B[i-1]);

//fxym[0] = fxym[0] / 2;

for(i = 1; i <= n; i++)

fxym[i] = (fxym[i] - a[i]*fxym[i-1]) / (2 - a[i]*B[i-1]); for(i = n-1; i >= 0; i--)

fxym[i] = fxym[i] - B[i]*fxym[i+1];

}

void printout(int n);

int main()

{

int n,i; char ch;

do

三次样条插值 C C++ 程序

{

cout<<"请输入已知断点个数:";

cin>>n;

for(i = 0; i <= n; i++)

{

cout<<"Please put in X"<<i<<':';

cin>>x[i]; //cout<<endl;

cout<<"Please put in Y"<<i<<':';

cin>>y[i]; //cout<<endl;

}

for(i = 0; i < n; i++) //求步长；其数组值较之x，y个数少一 h[i] = x[i+1] - x[i];

cout<<"Please 输入边界条件\n 1: 已知两端的一阶导数\n 2:两端的二阶导数已知\n 默认:自然边界条件\n";

int t;

float f0, f1;

cin>>t;

switch(t)

{

case 1:cout<<"Please put in Y0\' Y"<<n<<"\'\n";//显示数据为Y0'至Yn'，即断点的一阶导数

cin>>f0>>f1;

c[0] = 1; a[n] = 1;

fxym[0] = 6*((y[1] - y[0]) / (x[1] - x[0]) - f0) / h[0]; fxym[n] = 6*(f1 - (y[n] - y[n-1]) / (x[n] - x[n-1])) / h[n-1]; break;

case 2:cout<<"Please put in Y0\" Y"<<n<<"\"\n";//显示数据为Y0"至Yn"，即断点的二阶导数

cin>>f0>>f1;

c[0] = a[n] = 0;

fxym[0] = 2*f0; fxym[n] = 2*f1;

break;

default:cout<<"不可用\n";//待定

};//switch

for(i = 1; i < n; i++)

fxym[i] = 6 * f(i-1, i, i+1);//调用差分函数(only!) for(i = 1; i < n; i++)

{

a[i] = h[i-1] / (h[i] + h[i-1]);

c[i] = 1 - a[i];

}

a[n] = h[n-1] / (h[n-1] + h[n]);

cal_m(n);//调用弯矩函数（only！）

三次样条插值 C C++ 程序

cout<<"\n输出三次样条插值函数：\n";

printout(n);//调用求解三次样条插值函数；函数输出

cout<<"Do you to have anther try ? y/n :";

cin>>ch;

}

while(ch == 'y' || ch == 'Y');

return 0;

}

void printout(int n)/***************求三次样条插值函数（因已知断点个数而异）***********************/

{

cout<<setprecision(6);//通过操作器setprecision()设置有效位数；其为头文件<iomanip.h>所包含；括号内为参数。

for(int i = 0; i < n; i++)//所输出函数个数由所设断点个数而定 {

cout<<i+1<<": ["<<x[i]<<" , "<<x[i+1]<<"]\n"<<"\t";

float t = fxym[i]/(6*h[i]);

if(t > 0)

cout<<t<<"*("<<x[i+1]<<" - x)^3";

else

cout<<-t<<"*(x - "<<x[i+1]<<")^3";

t = …… 此处隐藏：5442字，全部文档内容请下载后查看。喜欢就下载吧 ……