【C/C++实现连分式插】

最新推荐文章于 2023-02-25 10:36:52 发布

Yasen.M

最新推荐文章于 2023-02-25 10:36:52 发布

阅读量413

点赞数

分类专栏：常用算法（C++描述）文章标签： c++ c语言

本文链接：https://blog.csdn.net/qq_45449625/article/details/127002039

版权

常用算法（C++描述）专栏收录该内容

13 篇文章 0 订阅

订阅专栏

【功能】
给定n个结点 $x_{i}$ (i=0,1,…,n-1)上的函数值 $y_{i}=f(x_{i})$ ,用连分式插值法计算指定值点t处的函数近似值z=f(t)。
【方法说明】
设给定的结点为 $x_{i}<x_{2}...<x_{n-1}<x_{n}$ ,其相应的函数值为 $y_{1},y_{2},...,y_{n-1},y_{n}$ ,则可以构造一个n-1节连分式

$\varphi (x)=b_{0}+\frac{x-x{0}}{b_{1}+\frac{x-x_{1}}{b_{2}+\frac{...}{...+\frac{x-x_{n-2}}{b_{n-1}}}}}$

计算 $b_{j}=\varphi _{j}(x_{i})$ (j=0,1,…,n-1)的递推公式如下:

$b_{0}=\varphi_{0}\left(x_{0}\right)=f\left(x_{0}\right) \\ \left\{\begin{array}{l} \varphi_{0}\left(x_{j}\right)=f\left(x_{j}\right) \\ \varphi_{k+1}\left(x_{j}\right)=\frac{x_{j}-x_{k}}{\varphi_{k}\left(x_{j}\right)-b_{k}}, \quad k=0,1, \cdots, j-1 \\ b_{j}=\varphi_{j}\left(x_{j}\right) \end{array}\right.$
在实际进行递推计算时,考虑到各中间的 $\varphi _{k}(x_{i})$ 值只被下一次使用。因此,可以将上述递推公式改写成如下形式:

$\left\{\begin{array}{l} u=f\left(x_{j}\right) \\ u=\frac{x_{j}-x_{k}}{u-b_{k}}, \quad k=0,1, \cdots, j-1 \\ b_{j}=u \end{array}\right.$
其中 $b_{0}=f(x_{0})$ 。
在实际进行插值计算时，一般在指定插值点t的前后各取4个结点就够了。此时，计算7节连分式的值c(t)作为插值点t处的函数近似值。

//【函数程序】
//连分式逐步插值.cpp
#include <cmath>	
#include <iostream> 
using namespace std;
//计算函数连分式新一节
//x 存放结点值x[0]~x[j]
//y 存放结点函数值y[0]~y[j]
//b存放连分式中的参数b[0]~b[j-1],返回时新增加b[j]
//j 连分式增加的节号。即本函数计算新增加的b[j]
void funpqi (double x[],double y[l,double b[l,int j)
{
	int k,flag=0; 
	double u;
	u=y[j];
	for(k=0;(k<j)&&(flag==0);k++)
		{
			if((u-b[k])+1.0-=1.0)flag-l; else
			u=(x[j]-x[k])/(u-b[k]);
		}
	if(flag==1)	u=1.0e+35; 
	b[j]=u;
	return;
}
//计算函数连分式值
//x 存放n个结点值x[0]~x[n-1]
//b 存放连分式中的n+1个参数b[0]~b[n]//n 
//n 连分式的节数(注意:常数项b[0]为第0节)
//t 自变量值
//程序返回t处的函数连分式值 
double funpqv(double x[],double b[],int n,double t)
{
	int k; double u; u=b[n];
		for(k=n-1;k>=0;k--)
		{
			if(fabs(u)+1.0==1.0)
				u=1.0e+35*(t-x[k])/fabs(t-x[k]); 
			else
				u=b[k]+(t-x[k])/u;
		}
	return(u);
}
//连分式逐步插值
//x[n]	存放结点值x[0]~x[n-1]	
//y[n]	存放结点函数值y[0]~y[n-1]	
//n	数据点个数。实际插值时最多取离插值点t最近的8个点	
//eps	控制精度要求	
//t	插值点值	
//返回插值点t处的连分式函数值
double funpq(double x[],double y[l,int n,double eps,double t)
{
	int i,j,k,l,m; 
	double p,q,u;
	//最多取离插值点t最近的8个点
	cout<<"最多取插值点最近的K个点"; 
	double b[8],xx[8],yy[8]; 
	p=0.0;
	if(n<1) return(p);	//结点个数不对，返回
	if(n==1)	{p=y[0]; return(p);}	//只有一个结点,取值返回	
	m=8;
	//最多取 8个点
	if (m>n)m=n;
	if(t<=x[0]) k=1;
	//第一个结点离插值点最近
	else if (t>=x[n-1]) k=n;
	//最后一个结点离插值点最近
	else
	{
		k=1;j=n;
		while((k-j!=1)&&(k-j!=-1))//二分法寻找离插值点最近的点
		{
			l=(k+j)/2;
			if(t<x[l-1])j=l; 
			else k=l;
		}
		if(fabs(t-x[1-1])>fabs(t-x[j-1])) k=j;
	}
	j=1; l=0;
	for (i=1;i<=m;i++)	//从数据表中取m个结点	
	{
		k=k+j*l;
		if((k<1)||(k>n))
		{
			l=l+1;j=-j;k=k+j*l;
		}
		xx[i-1]=x[k-1];
		yy[i-1]=y[k-1];
		l=l+1;j=-j;
	}
	j=0;
	b[0]=yy[0];
	p=b[0]; 
	u=1.0+eps;
	while((j<m-1)&&(u>=eps))
	{
		j=j+1;
		funpqj(xx,yy,b,j);	//计算新一节的b[j]	
		q=funpqv(xx,b,j,t);	//计算函数连分式在t处的函数值	
		u=(fabs(q-p)); p=q;
	}
	return(p);
}