分段线性插值.doc
插值法——线性分段插值

插值法——线性分段插值 1.插值函数%%分段线性插值function PLI = Piecewise_linear_interpolation(X,f,precision)[m,n] = size(X);a = min(X);b = max(X);X = sort(X);F = subs(f,X);for k = 1:n-1B = Basic_fun(X,k);I = B(1)*F(k)+B(2)*F(k+1);PLI{1,k} = [X(k),X(k+1)];PLI{2,k} = I;t{k} = X(k):(X(k+1)-X(k))/precision:X(k+1);T{k} = subs(I,t{k});Y_real{k} = subs(f,t{k});endfor k = 1:n-1t_((precision+1)*(k-1)+1:(precision+1)*k) = t{k};T_((precision+1)*(k-1)+1:(precision+1)*k) = T{k};Y_real_((precision+1)*(k-1)+1:(precision+1)*k)= Y_real{k};endh = figure;set(h,'color','w');plot(X,F,'r*',t_,T_,'g',t_,Y_real_,'b');xlabel('x shaft');ylabel('y shaft');legend('F:节点对应函数值','T:分段线性插值函数图像','Y_real:真实函数图像');title('分段线性插值');grid onend 2.基函数%%基函数,max(X)>k>0function BF = Basic_fun(X,k)X = sort(X);syms x;BF(1) = (x-X(k+1))/(X(k)-X(k+1));BF(2) = (x-X(k))/(X(k+1)-X(k));end 3.拟合值函数%%线性插值拟合值function LIV = Linear_interpolation_value(X,f,precision,x_value)[m,n] = size(X);a = min(X);b = max(X);X = sort(X);Answer = Piecewise_linear_interpolation(X,f,precision);for i = 1:n-1if x_value >= X(i) && x_value <= X(i+1)s = i;endendLIV{1,1} = '线性插值拟合值';LIV{2,1} = vpa(subs(Answer{2,s},x_value),6);LIV{1,2} = '真实值';LIV{2,2} = vpa(subs(f,x_value),6);LIV{1,3} = '误差';LIV{2,3} = abs(LIV{2,1}-LIV{2,2});end 4.例⼦clear allclcX = -5:1:5;syms x;f = - 0.08858*x^8 + 3.694*x^7 - 64.7*x^6 + 617.8*x^5 - 3490.0*x^4 + 11820.0*x^3 - 23150.0*x^2 + 23580.0*x - 9319.0; precision = 200;%%分段线性插值disp('分段线性插值');Piecewise_linear_interpolation(X,f,precision) 结果分段线性插值S =2×10 cell 数组列 1 ⾄ 4{1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×1 sym } {1×1 sym } {1×1 sym } {1×1 sym }列 5 ⾄ 8{1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×1 sym } {1×1 sym } {1×1 sym } {1×1 sym }列 9 ⾄ 10{1×2 double} {1×2 double}{1×1 sym } {1×1 sym }>> S{2,:}ans =(227077586881*x)/50000 + 37695704689/2500ans =(3983468847*x)/2000 + 60987657739/12500ans =(7723057429*x)/10000 + 30518164433/25000ans =(2518396259*x)/10000 + 4494858583/25000ans =(3136314129*x)/50000 - 9319ans =(465835271*x)/50000 - 9319ans =(422501*x)/10000 - 1113617/25000ans =4111433/25000 - (622509*x)/10000ans =- (271*x)/80 - 151661/12500ans =2072089/2500 - (10681481*x)/50000 图像如下。
分段线性插值

分段线性插值分段线性插值是一种在机器学习、数学、信号处理等领域中广泛应用的方法。
分段线性插值的主要目的是为漏洞、持续时间等数据展示提供更好的视觉效果,同时也可以使数据更容易进行处理。
在分段线性插值中,每一段数据都可以看作是一条直线段。
通过在相邻数据点之间插入一条直线来实现插值。
每个数据点或任意数段可以称为一个插值区间,插值区间内部的数据点都采用一条直线进行插值,直线的斜率由插值区间上下数据点构成。
例如:在一个区间(x1,y1)和(x2,y2)之间进行插值,其中x1<x<x2。
那么,我们可以使用线性公式y = mx + b来估计数据点的y值。
方程中m是插值区间的斜率,通过公式m = (y2-y1)/(x2-x1)计算。
而b是在插值区间x1和x2之间的截距,通过公式b = y1 - m x1计算。
最后,我们就可以通过已知的数据点,估计同一段中任意点的y值。
下面我们通过一个实例来进一步解释分段线性插值的应用。
比如我们有一组工作时间数据如下:|年份| 工作时间 ||----|----|| 2010 | 6.5 || 2011 | 7.0 || 2013 | 7.5 || 2015 | 8.0 |目前,我们需要在2012年估计工作时间。
首先,我们需要找到分段线性插值的区间。
2012年的数据点在2011年和2013年之间。
因此,我们可以使用2011年和2013年之间的数据点进行插值。
然后,通过计算斜率来确定m和b的值。
斜率可以通过公式m = (y2-y1)/(x2-x1)来计算。
2011年和2013年的工作时间分别是7.0和7.5,年份分别是2011和2013。
因此,斜率为:(7.5-7.0)/(2013年-2011年)= 0.25/2 = 0.125插值区间的y截距b可以通过公式b = y1 - m x1来计算。
这使得我们可以计算出截距:接下来,我们就可以使用斜率和截距来计算2012年的工作时间,这将是我们所需的数据点的估计值:y = mx + b= 0.125 * 2012 + 258.375= 259.875。
计算方法 1.3 分段线性插值

x x x x i 1 i ˆ S ( x ) y y , x x x 1 i i 1 i i 1 x x x x i i 1 i 1 i 于是, S1 ( x ) 是在 [ a , b ] 上是连续函数。
x [xj , xj ] 1 x [xj, xj 1] 其他
2)在插值节点 x 0 上,插值基为:
2 ( x x ) l ( x ) x [ x ,x ] 0 0 0 1 B ( x ) 0 0 其他 3)在插值节点 x n 上,插值基为:
2 ( x x ) l ( x ) n n B ( x ) n 0
1
左,右连接起来!
x
j1
xj
x
j1
2 2 H ( x ) 1 2 l ( x ) l ( x ) y 1 2 l ( x ) l ( x ) y 3 j 1 j j j j 1 j 1 2 2 ( x x ) l ( x ) y ( x x ) l ( x ) y j j j j 1j 1 j 1
k axb
提示:类似于前面的误差估计。 几点说明:
1)只要节点间距充分小,插值法总能获得所要求的精度。 2)局部性。如果修改某个数据,则插值曲线仅在某个局部范围内受影响。
插值节点 x 上,取值为 0 .即 k,k j 1 lj (x k ) 0 k j k j
2 )在每个小区间 [x 上,插值基 lj (x )都是线性函数 . i, x i 1]
基于以上两方面,我们观察
1
右 左
x
j1
线性分段插值函数

数值分析实验报告任课教师:马季骕班级:11级计算机科学与技术1实验目的及要求2程序的源代码3实验操作4实验结果及分析1实验目的及要求学会使用分段线性插值的方法求原函数的逼近函数。
随着插值节点的增加,插值多项式的插值多项式的次数也增加,而对于高次的插值容易带来剧烈的震荡,带来数值的不稳定(Runge现象)。
为了既要增加插值的节点,减小插值的区间,以便更好的逼近插值函数,又要不增加插值多项式的次数以减少误差,可采用分段线性插值.当给出了n+1个节点上f(x)的一张函数表后,用分段线性插值法求一个函数φ(x),并满足φ(x)是一个不超过n次的多项式;当插值节点取的足够多时逼近函数φ(x)能够很好的逼近被逼近函数f(x)。
而插值函数φ(x)的次数就会相应地升高,高次的插值多项式收敛到相应的被逼近函数。
所取的n值越大,逼近的越准确。
2程序源代码3// 数值分析Dlg.cpp : implementation file4//56#include "stdafx.h"7#include "数值分析.h"8#include "数值分析Dlg.h"910#ifdef _DEBUG11#define new DEBUG_NEW12#undef THIS_FILE13static char THIS_FILE[] = __FILE__;14#endif1516/////////////////////////////////////////////////////////////////////////// //17// CAboutDlg dialog used for App About1819class CAboutDlg : public CDialog20{21public:22CAboutDlg();2324// Dialog Data25//{{AFX_DATA(CAboutDlg)26enum { IDD = IDD_ABOUTBOX };27//}}AFX_DATA2829// ClassWizard generated virtual function overrides30//{{AFX_VIRTUAL(CAboutDlg)31protected:32virtual void DoDataExchange(CDataExchange* pDX); // DDX/DDV support33//}}AFX_VIRTUAL3435// Implementation36protected:37//{{AFX_MSG(CAboutDlg)38//}}AFX_MSG39DECLARE_MESSAGE_MAP()40};4142CAboutDlg::CAboutDlg() : CDialog(CAboutDlg::IDD)43{44//{{AFX_DATA_INIT(CAboutDlg)45//}}AFX_DATA_INIT46}4748void CAboutDlg::DoDataExchange(CDataExchange* pDX)49{50CDialog::DoDataExchange(pDX);51//{{AFX_DATA_MAP(CAboutDlg)52//}}AFX_DATA_MAP53}5455BEGIN_MESSAGE_MAP(CAboutDlg, CDialog)56//{{AFX_MSG_MAP(CAboutDlg)57// No message handlers58//}}AFX_MSG_MAP59END_MESSAGE_MAP()6061/////////////////////////////////////////////////////////////////////////////62// CMyDlg dialog6364CMyDlg::CMyDlg(CWnd* pParent /*=NULL*/)65: CDialog(CMyDlg::IDD, pParent)66{67//{{AFX_DATA_INIT(CMyDlg)68// NOTE: the ClassWizard will add member initialization here69//}}AFX_DATA_INIT70// Note that LoadIcon does not require a subsequent DestroyIcon in Win32 71m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME);72}7374void CMyDlg::DoDataExchange(CDataExchange* pDX)75{76CDialog::DoDataExchange(pDX);77//{{AFX_DATA_MAP(CMyDlg)78// NOTE: the ClassWizard will add DDX and DDV calls here79//}}AFX_DATA_MAP80}8182BEGIN_MESSAGE_MAP(CMyDlg, CDialog)83//{{AFX_MSG_MAP(CMyDlg)84ON_WM_SYSCOMMAND()85ON_WM_PAINT()86ON_WM_QUERYDRAGICON()87ON_BN_CLICKED(IDC_NEW, OnNew)88ON_BN_CLICKED(IDC_LAG, OnLag)89ON_BN_CLICKED(IDC_LINE, OnLine)90ON_BN_CLICKED(IDC_HER, OnHer)91//}}AFX_MSG_MAP92END_MESSAGE_MAP()9394/////////////////////////////////////////////////////////////////////////// //95// CMyDlg message handlers9697BOOL CMyDlg::OnInitDialog()98{99CDialog::OnInitDialog();100101// Add "About..." menu item to system menu.102103// IDM_ABOUTBOX must be in the system command range.104ASSERT((IDM_ABOUTBOX & 0xFFF0) == IDM_ABOUTBOX);105ASSERT(IDM_ABOUTBOX < 0xF000);106107CMenu* pSysMenu = GetSystemMenu(FALSE);108if (pSysMenu != NULL)109{110CString strAboutMenu;111strAboutMenu.LoadString(IDS_ABOUTBOX);112if (!strAboutMenu.IsEmpty())113{114pSysMenu->AppendMenu(MF_SEPARATOR);115pSysMenu->AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu);116}117}118119// Set the icon for this dialog. The framework does this automatically120// when the application's main window is not a dialog121SetIcon(m_hIcon, TRUE); // Set big icon122SetIcon(m_hIcon, FALSE); // Set small icon123124// TODO: Add extra initialization here125126return TRUE; // return TRUE unless you set the focus to a control127}128129void CMyDlg::OnSysCommand(UINT nID, LPARAM lParam)130{131if ((nID & 0xFFF0) == IDM_ABOUTBOX)132{133CAboutDlg dlgAbout;134dlgAbout.DoModal();135}136else137{138CDialog::OnSysCommand(nID, lParam);139}140}141142// If you add a minimize button to your dialog, you will need the code below 143// to draw the icon. For MFC applications using the document/view model, 144// this is automatically done for you by the framework.145146void CMyDlg::OnPaint()147{148if (IsIconic())149{150CPaintDC dc(this); // device context for painting151152SendMessage(WM_ICONERASEBKGND, (WPARAM) dc.GetSafeHdc(), 0);153154// Center icon in client rectangle155int cxIcon = GetSystemMetrics(SM_CXICON);156int cyIcon = GetSystemMetrics(SM_CYICON);157CRect rect;158GetClientRect(&rect);159int x = (rect.Width() - cxIcon + 1) / 2;160int y = (rect.Height() - cyIcon + 1) / 2;161162// Draw the icon163dc.DrawIcon(x, y, m_hIcon);164}165else166{167CDialog::OnPaint();168}169}170171// The system calls this to obtain the cursor to display while the user drags 172// the minimized window.173HCURSOR CMyDlg::OnQueryDragIcon()174{175return (HCURSOR) m_hIcon;176}177178void CMyDlg::OnNew()179{180// TODO: Add your control notification handler code here181 int x00=225,y00=225,i,j;182double x;183184CDC *pDC=GetDC();185pDC->SetMapMode(MM_LOMETRIC);186187pDC->SetViewportOrg(x00,y00);188189//画坐标轴与原函数190for(i=-650; i<=650; i++)191{192pDC->SetPixel(i,0,RGB(0,0,0));193pDC->SetPixel(0,i,RGB(0,0,0));194}195196for(x=-1; x<=1; x+=0.001)197{198double j=400.0/(1+25*x*x);199pDC->SetPixel(x*500,j,RGB(255,0,0));200}201202}203204void CMyDlg::OnLag()205{206// TODO: Add your control notification handler code here 207int x00=225,y00=225,i,j;208double x;209210CDC *pDC=GetDC();211pDC->SetMapMode(MM_LOMETRIC);212213pDC->SetViewportOrg(x00,y00);214215//画坐标轴216for(i=-650; i<=650; i++)217{218pDC->SetPixel(i,0,RGB(0,0,0));219pDC->SetPixel(0,i,RGB(0,0,0));220}221222223double yx[]={-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1}; 224225// 拉格朗日差值的函数226double yy[12],lx[12],ly[12];227double l_fenzi[12],l_fenmu[12];228 double l_x,l_y;229for(i=0; i<=10; i++)230{231yy[i]=1.0/(1+25*yx[i]*yx[i]);232}233234for(i=0; i<=10; i++)235{236l_fenmu[i]=1.0;237for(j=0; j<=10; j++)238{239if(i!=j)240l_fenmu[i]=l_fenmu[i]*(yx[i]-yx[j]);241}242}243 double qq,pp;244for(qq=-1; qq<=1; qq+=0.0003)245{246for(i=0; i<=10; i++)247{248l_fenzi[i]=1.0;249for(j=0; j<=10; j++)250{251if(i!=j)252 l_fenzi[i]=l_fenzi[i]*(qq-yx[j]);253}254}255pp=0;256for(i=0; i<=11; i++)257{258 pp=pp+1.0/(1+25*yx[i]*yx[i])*l_fenzi[i]/l_fenmu[i]; 259}260pDC->SetPixel(qq*500,pp*390+5,RGB(150,255,0));261}262263}264265void CMyDlg::OnLine()266{267// TODO: Add your control notification handler code here 268int x00=225,y00=225,i,j;269double x;270271CDC *pDC=GetDC();272pDC->SetMapMode(MM_LOMETRIC);273274pDC->SetViewportOrg(x00,y00);275276//画坐标轴与原函数277for(i=-650; i<=650; i++)278{279pDC->SetPixel(i,0,RGB(0,0,0));280pDC->SetPixel(0,i,RGB(0,0,0));281}282283double yx[]={-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1}; 284double yy[14];285for(i=0; i<=10; i++)286{287yy[i]=1.0/(1+25*yx[i]*yx[i]);288}289290// 线性分段差值的图像291CPen pen;292293CPen*oldpen;294295pen.CreatePen(PS_SOLID,5,RGB(0,0,0));296297oldpen=pDC->SelectObject(&pen);298for(i=0; i<10; i++)299{300 pDC->MoveTo(yx[i]*480,yy[i]*400);301pDC->LineTo(yx[i+1]*480,yy[i+1]*400);302 }303}304305306void CMyDlg::OnHer()307{308// TODO: Add your control notification handler code here 309int x00=225,y00=225,i,j;310double x;311312CDC *pDC=GetDC();313pDC->SetMapMode(MM_LOMETRIC);314315pDC->SetViewportOrg(x00,y00);316317//画坐标轴与原函数318for(i=-650; i<=650; i++)319{320pDC->SetPixel(i,0,RGB(0,0,0));321pDC->SetPixel(0,i,RGB(0,0,0));322}323324 double yx[]={-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1};325326 double yy[12];327 for(i=0; i<=10; i++)328 {329 yy[i]=1.0/(1+25*yx[i]*yx[i]);330 }331332 //分段三次Hermite 差值的函数333334 double x0,x1,yd1,yd0,y1,y0;335 for(i=0; i<10; i++)336 {337338 x0=yx[i],x1=yx[i+1];339 y0=1.0/(1+25*x0*x0);340 y1=1.0/(1+25*x1*x1);341 yd0=-(50*x0)*1.0/((1+25*x0*x0)*(1+25*x0*x0));342 yd1=-(50*x1)*1.0/((1+25*x1*x1)*(1+25*x1*x1));343344 for(double qq=x0; qq<x1; qq+=0.005)345 {346 double pp= y0*(1+2*(qq-x0)/(x1-x0)) * (qq-x1)/(x0-x1) *(qq-x1)/(x0-x1)347 +y1*(1+2*(qq-x1)/(x0-x1)) * (qq-x0)/(x1-x0) * (qq-x0)/(x1-x0) 348 +yd0*(qq-x0) * (qq-x1)/(x0-x1) * (qq-x1)/(x0-x1)349 +yd1*(qq-x1) * (qq-x0)/(x1-x0) * (qq-x0)/(x1-x0);350 pDC->SetPixel(qq*500,pp*400,RGB(0,255,255));351 }352 }353354 }程序操作算法分析:1. 分段线性插值的算法思想:分段线性插值需要在每个插值节点上构造分段线性插值基函数)(x l j ,然后再作它们的线性组合。
分段线性插值

三次样条插值
插值函数S(x): (1)S(x)在字区间[xi,xi+1]的表达式Si(x)都是次数不高 于3的多项式; (2) Si(x)=yi; (3) S(x)在整个区间[a,b]上有连续的二阶导数。
hk xk 1 xk
yk 1 yk hk S ( x) yk 1 [ (2mk 1 mk )](x xk 1 ) hk 6
[fnodes,minnq,rnw,rnq,ifail]=e01sef(x,y,z) [cz(i,j),ifail]=e01sff(x,y,z,rnw,fnodes,cx(I),cy(j))
数据拟合
已知平面上n个点,寻求函数f(x)在某种准则下与数据点 最为接近,即曲线拟合得好。
‘nearest’-----最近邻点插值;
‘linear’-------线性插值;-----缺省 ‘spline’-------三次样条插值; ‘cubic’--------三次插值。
散乱数据插iddata(x,y,z,cx,cy,’method’)
返回在网格(cx,cy)处的函数值,注意此处cx行向 量,cy列向量。 (2)插值函数e01sef和e01sff 基于Shephard插值法,两个函数必须同时使用。
mk S ( xk ) mk 1 mk mk 1 2 hk ( x xk 1 ) ( x xk 1 )3 k 2 6hk 1 hk 1 hk y y y yk 1 6 k ( k 1 k k ) n+1个未知量,n-1个方程组 hk 1 hk hk hk 1
二维插值
通过全部已知节点f(xi,yi)=zij(i=0,1,2, …m;j=1,1,2, …n), 构造一二元函数z=f(x,y),然后再利用f(x,y)插值,即 z*= f(x*,y *);分为网格节点插值和散乱节点插值两种。
分段线性插值

摘要用函数来表示变量间的数量关系广泛应用于各学科领域,但是在实际问题中,往往是通过实验、观测以及计算等方法,得到的是函数在一些点上的函数值。
如何通过这些离散数据找到函数的一个满足精度要求且便于使用的近似表达式,是经常遇到的问题。
对于这类问题我们解决的方法为插值法,而最常用也最简单的插值方法就是多项式插值。
当然用插值法得到的近似表达式必须满足插值条件即假设给定了n+1个点的自变量的值以及函数值,近似函数必须要过这n+1(x)通个点。
多项式插值,从几何角度看,就是寻求n次代数曲线y=Pn过n+1个点作为f(x)的近似。
但是随着插值节点个数的增加,高次插值多项式的近似效果并不理想。
根据大量实验得出,在进行高次多项式插值时,会出现龙格现象。
龙格(Runge)现象即当n趋于无穷大时,x在某一邻域内,f(x)收敛,而在这个区域外f(x)发散。
因此,为了解决这样的一个问题,我们可以通过缩小插值区间的办法达到减小误差的目的,所以本实验将针对低次分段插值多项式来做具体的讨论和学习。
关键词:龙格现象分段差值1、实验目的1)通过对分段线性插值算法程序的编写,提高自己编写程序的能力2)体会分段线性插值是如何消除龙格现象的。
3)用实验报告的形式展现,提高自己在写论文方面的能力2、算法理论设在节点处的函数值为,i=0,1,,n。
为了提高近似程度,可以考虑用分段线性插值来逼近原函数,这时的插值函数为分段函数:在区间上的线性函数为误差为:易见,是平面上以点为节点的折线,有如下的特点:1.在上为次数不超过一次的多项式;2.;3.;如果,由线性插值的误差公式得到令,则有关于整体误差:可以按如下方式考虑,若记则对任一都有于是,当时,说明分段线性插值收敛于。
3、数值算例已知点坐标如下表所示:xiyi用分段线性插值法,求解当x为时,对应y的值解:具体程序如下所示:#include ""float Fdline(float x[],float y[],float x1,int len){int i=0;float s=0;for(i=0;i<len-1;i++){if(x1>=x[i] && x1<x[i+1])break;}s=(x1-x[i])/(x[i-1]-x[i])*y[i-1]+(x1-x[i-1])/(x[i]-x[i-1])* y[i];return s;}float Fdline(float x[],float y[],float x1,int len);void main(){float x[]={,,,,};float y[]={,,,,};int len=sizeof(x)/sizeof(x[0]);float x1=0;float s=0;printf("请输入要求解的x1的值:\n");scanf("%f",&x1);s=Fdline(x,y,x1,len);printf("经过分段三次Hermite插值的结果为:\n");printf("%f\n",s);}运行结果:5、对结果进行分析根据分段线性插值的原理,可以看出分段线性插值虽然有很好的收敛性质,但却不是光滑的,所以线性插值的结果和实际的结果差距较大。
分段线性插值法

《数值分析》实验报告实验序号:实验五 实验名称: 分段线性插值法1、 实验目的:随着插值节点的增加,插值多项式的插值多项式的次数也增加,而对于高次的插值容易带来剧烈的震荡,带来数值的不稳定(Runge 现象)。
为了既要增加插值的节点,减小插值的区间,以便更好的逼近插值函数,又要不增加插值多项式的次数以减少误差,可采用分段线性插值。
2、 实验内容:求一个函数ϕ(x )用来近似函数f (x ),用分段线性插值的方法来求解近似函数ϕ(x )并画出近似函数图像及原函数图像。
设在区间[a,b]上,给定n+1个插值节点b x x x x a n =<<<<=...210与相应的函数值n y y y ,...,,10,求一个插值函数)(x ϕ,满足以下条件:(1)),...,2,1,0()(n j y x j j ==ϕ; (2) )(x ϕ在每一个小区间[1,+j j x x ]上就是线性函数。
对于给定函数11-,2511)(2≤≤+=x x x f 。
在区间[]11-,上画出f (x )与分段线性插值函数)(x ϕ的函数图像。
1. 分段线性插值的算法思想:分段线性插值需要在每个插值节点上构造分段线性插值基函数)(x l j ,然后再作它们的线性组合。
分段线性插值基函数的特点就是在对应的插值节点上函数值取 1,其它节点上函数值取0。
插值基函数如下:⎪⎩⎪⎨⎧≤≤--=其它 ,0,)(101010x x x x x x x x l ⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧≤<--≤≤--=+++---其它 ,0,,)(111111j j j j j j j j j j j x x x x x x x x x x x x x x x l⎪⎩⎪⎨⎧≤≤--=---其它 ,0,)(111n n n n n n x x x x x x x x l设在节点a≤x0<x1<…≤b=f(xi),(i=0,1,2,…,n)求折线函数L(x)满足:(1)L(x)∈C[a,b](2)L(x[i]=y[i])(3)L(x)在每个小区间(x[i],x[i+1])上就是线性插值函数¢(x)叫做区间[a,b]上对数据(x[j],y[j])(j=0,1,2,…,n)的分段区间函数。
分段线性插值

1、实验目的:设在区间[a,b]上,给定n+1个插值节点a=X0<X1<X2<。
<Xn=b和相应的函数值y0,y1,y2,。
,yn ,求一个插值函数φ(x),具有下面性质:(1)φ(x)=yj(j=0,1,2,…,n)(2)φ(x)在每一个小区间[xj,yj+1]上是线性函数。
2、算法分析:●分段线性插值的算法思想:分段线性插值需要在每个插值节点上构造分段线性插值基函数,然后再作它们的线性组合。
分段线性插值基函数的特点是在对应的插值节点上函数值取1,其它节点上函数值取0。
插值基函数如下:有了基函数就可以直接写出分段线性插值函数的表达式。
●具体程序设计:for(i=0;i<=20;i++) //选取节点{ax[i]=-1+((2/20.0)*i); //选取上的21个对称的节点ay[i]=1.0/(1+25*ax[i]*ax[i]);}x1=-1;while(x1<=1){m=0;x1=x1+0.00001;for(i=0;i<=20;i++){if(i==0&&x1>=-1&&x1<=-0.9){n=ay[0]*((x1-ax[1])/(ax[0]-ax[1])); //计算并和相应的函数值组合m=n+m;}else{if(x1>=ax[i-1]&&x1<=ax[i]) //计算并和相应的函数值组合{n=ay[i]*((x1-ax[i-1])/(ax[i]-ax[i-1]));m=n+m;}else{if(x1>ax[i]&&x1<=ax[i+1]){n=ay[i]*((x1-ax[i+1])/(ax[i]-ax[i+1]));m=m+n;}else{if(i==19&&x1>ax[19]&&x1<=ax[20]){n=ay[20]*((x1-ax[19])-(ax[20]-ax[19]));//计算并和相应的函数值组合m=n+m;}}}}}}3、实验结果截图:在[-1,1]区间上选取了21个等分节点的分段线性插值函数的图像如下:4、程序代码// SHIYAN456View.cpp : implementation of the CSHIYAN456View class //#include "stdafx.h"#include "SHIYAN456.h"#include "SHIYAN456Doc.h"#include "SHIYAN456View.h"#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif/////////////////////////////////////////////////////////////////////////////// CSHIYAN456ViewIMPLEMENT_DYNCREATE(CSHIYAN456View, CView)BEGIN_MESSAGE_MAP(CSHIYAN456View, CView)//{{AFX_MSG_MAP(CSHIYAN456View)ON_COMMAND(ID_FFunction, OnFFunction)ON_COMMAND(ID_Lagrange, OnLagrange)ON_COMMAND(ID_Subsection, OnSubsection)ON_COMMAND(ID_Hermite, OnHermite)//}}AFX_MSG_MAP// Standard printing commandsON_COMMAND(ID_FILE_PRINT, CView::OnFilePrint)ON_COMMAND(ID_FILE_PRINT_DIRECT, CView::OnFilePrint)ON_COMMAND(ID_FILE_PRINT_PREVIEW, CView::OnFilePrintPreview) END_MESSAGE_MAP()/////////////////////////////////////////////////////////////////////////////// CSHIYAN456View construction/destructionCSHIYAN456View::CSHIYAN456View(){// TODO: add construction code here}CSHIYAN456View::~CSHIYAN456View(){}BOOL CSHIYAN456View::PreCreateWindow(CREATESTRUCT& cs)// TODO: Modify the Window class or styles here by modifying // the CREATESTRUCT csreturn CView::PreCreateWindow(cs);}/////////////////////////////////////////////////////////////////////////////// CSHIYAN456View drawingvoid CSHIYAN456View::OnDraw(CDC* pDC){CSHIYAN456Doc* pDoc = GetDocument();ASSERT_VALID(pDoc);// TODO: add draw code for native data hereAfxGetMainWnd()->SetWindowText("实验四五六函数图像");for(int k=650;k<=690;k++){pDC->SetPixel(k,55,RGB(255,0,0));pDC->SetPixel(k,85,RGB(0,255,0));pDC->SetPixel(k,115,RGB(0,0,255));pDC->SetPixel(k,145,RGB(0,255,255));}pDC->TextOut(700,50,"原函数图像");pDC->TextOut(700,80,"Lagrange插值函数图像");pDC->TextOut(700,110,"分段线性插值函数图像");pDC->TextOut(700,140,"Hermite插值函数图像");for(int i=6;i<=600;i++)pDC->SetPixel(400,i,RGB(0,0,0));pDC->TextOut(395,4,"y");for(int j=100;j<=700;j++)pDC->SetPixel(j,500,RGB(0,0,0));pDC->TextOut(700,500,"x");//pDC->MoveTo(400,500);}/////////////////////////////////////////////////////////////////////////////// CSHIYAN456View printingBOOL CSHIYAN456View::OnPreparePrinting(CPrintInfo* pInfo){// default preparationreturn DoPreparePrinting(pInfo);}void CSHIYAN456View::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/) {// TODO: add extra initialization before printing}void CSHIYAN456View::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/) {// TODO: add cleanup after printing}/////////////////////////////////////////////////////////////////////////////// CSHIYAN456View diagnostics#ifdef _DEBUGvoid CSHIYAN456View::AssertValid() const{CView::AssertValid();}void CSHIYAN456View::Dump(CDumpContext& dc) const{CView::Dump(dc);}CSHIYAN456Doc* CSHIYAN456View::GetDocument() // non-debug version is inline {ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CSHIYAN456Doc)));return (CSHIYAN456Doc*)m_pDocument;}#endif //_DEBUG/////////////////////////////////////////////////////////////////////////////// CSHIYAN456View message handlersvoid CSHIYAN456View::OnFFunction(){// TODO: Add your command handler code hereCClientDC dr(this);COLORREF rgb=RGB(255,0,0);double x1,y1,x,y;x1=-1.0;y1=1/(1+25*x1*x1);x=x1*200+400;y=-y1*200+500;while(x1<=1){dr.MoveTo(int(x),int(y));x1=x1+0.00001;y1=1/(1+25*x1*x1);x=x1*200+400;y=-y1*200+500;dr.SetPixel(int(x),int(y),rgb);}}void CSHIYAN456View::OnLagrange(){// TODO: Add your command handler code hereCClientDC dr(this);COLORREF rgb=RGB(0,255,0);int i,j;double x1=0,y1=0,x=0,y=0,m=0,n=0,ax[100],ay[100];for(i=0;i<=10;i++){ax[i]=-1+((2/10.0)*i);ay[i]=1.0/(1+25*ax[i]*ax[i]);}x1=-1;y1=1/(1+25*x1*x1);x=x1*200+400;y=-y1*200+500;dr.MoveTo(int(x),int(y));while(x1<=1){m=0;for(i=0;i<=10;i++){n=1;for(j=0;j<=10;j++){if(i!=j)n=((x1-ax[j])/(ax[i]-ax[j]))*n;}m=ay[i]*n+m;}x=x1*200+400;y=-m*200+500;dr.SetPixel(int(x),int(y),rgb);x1=x1+0.00001;}}void CSHIYAN456View::OnSubsection(){// TODO: Add your command handler code hereCClientDC dr(this);COLORREF rgb=RGB(0,0,255);int i;double x1=0,y1=0,x=0,y=0,m=0,n=0,ax[100],ay[100];for(i=0;i<=20;i++){ax[i]=-1+((2/20.0)*i);ay[i]=1.0/(1+25*ax[i]*ax[i]);}x1=-1;y1=1/(1+25*x1*x1);x=x1*200+400;y=-y1*200+500;dr.MoveTo(int(x),int(y));while(x1<=1){m=0;x1=x1+0.00001;for(i=0;i<=20;i++){if(i==0&&x1>=-1&&x1<=-0.9){n=ay[0]*((x1-ax[1])/(ax[0]-ax[1]));m=n+m;}else{if(x1>=ax[i-1]&&x1<=ax[i]){n=ay[i]*((x1-ax[i-1])/(ax[i]-ax[i-1]));m=n+m;}else{if(x1>ax[i]&&x1<=ax[i+1]){n=ay[i]*((x1-ax[i+1])/(ax[i]-ax[i+1]));m=m+n;}else{if(i==19&&x1>ax[19]&&x1<=ax[20]){n=ay[20]*((x1-ax[19])-(ax[20]-ax[19]));m=n+m;}}}}}x=x1*200+400;y=-m*200+500;dr.SetPixel(int(x),int(y),rgb);}}void CSHIYAN456View::OnHermite(){// TODO: Add your command handler code hereCClientDC dr(this);COLORREF rgb=RGB(0,255,255);int i;double x1=0,y1=0,x=0,y=0,m=0,n=0,tt=0,tt1=0,h=0,ax[100],ay[100],a[100];for(i=0;i<=10;i++){ax[i]=-1+((2/10.0)*i);ay[i]=1.0/(1+25*ax[i]*ax[i]);a[i]=(-50*ax[i])/((1+25*ax[i]*ax[i])*(1+25*ax[i]*ax[i]));}x1=-1;y1=1/(1+25*x1*x1);x=x1*200+400;y=-y1*200+500;dr.MoveTo(int(x),int(y));while(x1<=1){m=0;x1=x1+0.00001;for(i=0;i<=10;i++){if(i==0&&x1>=ax[0]&&x1<=ax[1]){n=(1+2*((x1-ax[0])/(ax[1]-ax[0])))*((x1-ax[1])/(ax[0]-ax[1]))*((x1-ax[1])/(ax[0]-ax[1]));tt=ay[0]*n;h=(x1-ax[0])*((x1-ax[1])/(ax[0]-ax[1]))*((x1-ax[1])/(ax[0]-ax[1]));tt1=a[0]*h;m=tt+tt1+m;}else{if(x1>=ax[i-1]&&x1<=ax[i]){n=(1+2*((x1-ax[i])/(ax[i-1]-ax[i])))*((x1-ax[i-1])/(ax[i]-ax[i-1]))*((x1-ax[i-1])/(ax[i]-ax[i-1]));tt=ay[i]*n;h=(x1-ax[i])*((x1-ax[i-1])/(ax[i]-ax[i-1]))*((x1-ax[i-1])/(ax[i]-ax[i-1]));tt1=a[i]*h;m=tt+tt1+m;}else{if(x1>ax[i]&&x1<=ax[i+1]){n=(1+2*((x1-ax[i])/(ax[i+1]-ax[i])))*((x1-ax[i+1])/(ax[i]-ax[i+1]))*((x1-ax[i+1])/(ax[i]-ax[i+1]));tt=ay[i]*n;h=(x1-ax[i])*((x1-ax[i+1])/(ax[i]-ax[i+1]))*((x1-ax[i+1])/(ax[i]-ax[i+1]));tt1=a[i]*h;m=tt+tt1+m;}else{if(i==9&&x1>ax[9]&&x1<=ax[10]){n=(1+2*((x1-ax[10])/(ax[9]-ax[10])))*((x1-ax[9])/(ax[10]-ax[9]))*((x1-ax[9])/(ax[10]-ax[9]));tt=ay[i]*n;h=(x1-ax[10])*((x1-ax[9])/(ax[10]-ax[9]))*((x1-ax[9])/(ax[10]-ax[9]));tt1=a[i]*h;m=tt+tt1+m;}}}}}x=x1*200+400;y=-m*200+500;dr.SetPixel(int(x),int(y),rgb);}}5、总结体会分段线性插值的方法克服了Lagrange插值法当节点不断加密时,构造的插值多项式的次数不断升高,高次多项式虽然是连续的,但是不一定都收敛到相应的被插函数而产生Runge现象。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
摘要
用函数来表示变量间的数量关系广泛应用于各学科领域,但是在实际问题中,往往是通过实验、观测以及计算等方法,得到的是函数在一些点上的函数值。
如何通过这些离散数据找到函数的一个满足精度要求且便于使用的近似表达式,是经常遇到的问题。
对于这类问题我们解决的方法为插值法,而最常用也最简单的插值方法就是多项式插值。
当然用插值法得到的近似表达式必须满足插值条件即假设给定了n+1个点的自变量的值以及函数值,近似函数必须要过这n+1
(x)通个点。
多项式插值,从几何角度看,就是寻求n次代数曲线y=P
n
过n+1个点作为f(x)的近似。
但是随着插值节点个数的增加,高次插值多项式的近似效果并不理想。
根据大量实验得出,在进行高次多项式插值时,会出现龙格现象。
龙格(Runge)现象即当n趋于无穷大时,x在某一邻域内,f(x)收敛,而在这个区域外f(x)发散。
因此,为了解决这样的一个问题,我们可以通过缩小插值区间的办法达到减小误差的目的,所以本实验将针对低次分段插值多项式来做具体的讨论和学习。
关键词:龙格现象分段差值
1、实验目的
1)通过对分段线性插值算法程序的编写,提高自己编写程序的能力
2)体会分段线性插值是如何消除龙格现象的。
3)用实验报告的形式展现,提高自己在写论文方面的能力
2、算法理论
设在节点处的函数值为,i=0,1,,n。
为了提高近似程度,可以考虑用分段线性插值来逼近原函数,这时的插值函数为分段函数:
在区间上的线性函数为
误差为:
易见,是平面上以点为节点的折线,有如下的特点:
1.在上为次数不超过一次的多项式;
2.;
3.;
如果,由线性插值的误差公式得到
令,则有
关于整体误差:
可以按如下方式考虑,若记则对任一都有
于是,当时,说明分段线性插值收敛于。
3、数值算例
0.5 0.7 0.9 1.1 1.3
x
i
y
0.4579 0.644 0.783 0.891 0.964
i
用分段线性插值法,求解当x为0.8时,对应y的值
解:具体程序如下所示:
#include "stdafx.h"
float Fdline(float x[],float y[],float x1,int len)
{
int i=0;
float s=0;
for(i=0;i<len-1;i++)
{
if(x1>=x[i] && x1<x[i+1])
break;
}
s=(x1-x[i])/(x[i-1]-x[i])*y[i-1]+(x1-x[i-1])/(x[i]-x[i-1])* y[i];
return s;
}
float Fdline(float x[],float y[],float x1,int len);
void main()
{
float x[]={0.5,0.7,0.9,1.1,1.3};
float y[]={0.479,0.644,0.783,0.891,0.964};
int len=sizeof(x)/sizeof(x[0]);
float x1=0;
float s=0;
printf("请输入要求解的x1的值:\n");
scanf("%f",&x1);
s=Fdline(x,y,x1,len);
printf("经过分段三次Hermite插值的结果为:\n");
printf("%f\n",s);
}
运行结果:
5、对结果进行分析
根据分段线性插值的原理,可以看出分段线性插值虽然有很好的收敛性质,但却不是光滑的,所以线性插值的结果和实际的结果差距较大。
通过用编程实现对上例的求解,可以看出结果较为准确,但是由于在计算机上计算,会存在计算误差。
6、参考文献
[1] 秦新强.数值逼近.西安:西安理工大学出版社,2010。