用vb实现利用三次样条插值函数进行编程

用vb实现利用三次样条插值函数进行编程

网友 cz5360 于提问

vb三次样条插值函数绘图

Dim X(1000) As Single, Y(1000) As Single

Dim u1(0 To 80000) As Single, v1(0 To 80000) As Single

Dim num As Long

Dim t As Integer

Private Declare Sub Sleep Lib "kernel32" (ByVal dwMilliseconds As Long)

Dim de As Integer

Dim ToInit As Boolean

Dim DownX As Single, DownY As Single

Sub Drawposi(Index As Integer)

Me.Picture1.ForeColor = 0

Me.Picture1.Line (0, Y(Index))-(1024, Y(Index))

Me.Picture1.Line (X(Index), 0)-(X(Index), 770)

End Sub

Function hypot(ByVal X As Single, ByVal Y As Single)

hypot = Sqr(X ^ 2 + Y ^ 2)

End Function

Sub MovePic(Index As Integer)

Dim i As Integer

X(Index) = Picture2(Index).Left + 4

Y(Index) = Picture2(Index).Top + 4

lblX.Caption = "X: " + CStr(CInt(X(Index)))

lblY.Caption = "Y: " + CStr(CInt(Y(Index)))

lblX.Refresh

lblY.Refresh

Me.Picture1.Cls

Me.Picture1.ForeColor = QBColor(10)

For i = 0 To t - 1

Me.Picture1.CurrentX = X(i) + 4

Me.Picture1.CurrentY = Y(i) + 4

Me.Picture1.Print i

Next i

End Sub

Private Sub Command1_Click()

Dim i As Long

'Picture1.Scale (0, 0)-(640, 550)

DrawWidth = 3

Picture1.Cls

'If Check1.Value Then Command2_Click

'X(0) = 1

'Y(0) = 1

'X(t - 1) = 638

'Y(t - 1) = 548

Picture1.ForeColor = QBColor(10)

For i = 0 To t - 1

Picture1.Line (X(i) - 1, Y(i) - 1)-(X(i) + 1, Y(i) + 1), QBColor(10), B Picture1.Print i

Next i

Picture1.ForeColor = QBColor(12)

DrawWidth = 1

tspLine t - 1, 2, 0, 0, 0, 0

Picture1.PSet (u1(0), v1(0))

For i = 1 To num - 1

Picture1.Line -(u1(i), v1(i))

'For de = 1 To 12000: Next de 'Sleep 1

Next i

Picture1.ForeColor = QBColor(10)

For i = 0 To t - 1

Picture1.Line (X(i) - 1, Y(i) - 1)-(X(i) + 1, Y(i) + 1), QBColor(10), B Picture1.Print i

Next i

End Sub

Private Sub Command2_Click()

Dim i As Integer

Randomize Timer

ToInit = Not ToInit

If ToInit Then

https://www.wendangku.net/doc/7d2289562.html,mand1.Enabled = False

https://www.wendangku.net/doc/7d2289562.html,mand2.Caption = "结束初始化"

Me.Cls

For i = 1 To t - 1

Load Me.Picture2(i)

Next i

For i = 0 To t - 1

Picture2(i).Left = X(i) - 4

Picture2(i).Top = Y(i) - 4

Picture2(i).Visible = True

Next i

Picture1.Cls

Me.Picture1.ForeColor = QBColor(10)

For i = 0 To t - 1

Picture1.Line (X(i) - 1, Y(i) - 1)-(X(i) + 1, Y(i) + 1), QBColor(10), B

Picture1.Print i

Next i

Else

https://www.wendangku.net/doc/7d2289562.html,mand1.Enabled = True

https://www.wendangku.net/doc/7d2289562.html,mand2.Caption = "开始初始化"

For i = 1 To t - 1

Unload Me.Picture2(i)

Next i

Me.Picture2(0).Visible = False

End If

Exit Sub

For i = 0 To t

X(i) = Rnd(1) * 500 + Rnd(1) * 50 + 12

Y(i) = Rnd(1) * 400 + Rnd(1) * 100 + 12

'X(i) = i * 20 + Rnd(1) * 10 + 12

'Y(i) = i * 10 + Rnd(1) * 300 + 22 - Rnd(1) * 200

Next i

End Sub

Sub tspLine(ByVal n As Integer, ByVal ch As Integer, ByVal tx1 As Single, ByVal tx2 As Single, ByVal ty1 As Single, ByVal ty2 As Single)

Dim a(1000) As Single, b(1000) As Single, c(1000) As Single, dX(1000) As Single, dY(1000) As Single

Dim qx(1000) As Single, qy(1000) As Single

Dim tt As Single, bx3 As Single, bx4 As Single, by3 As Single, by4 As Single Dim cx As Single, cy As Single, t(1000) As Single, px(1000) As Single,

py(1000) As Single

Dim u(3000) As Single, v(3000) As Single, i As Integer

num = 0

For i = 1 To n

t(i) = hypot(X(i) - X(i - 1), Y(i) - Y(i - 1))

Next i

Select Case ch

Case 0 '抛物条件

u(0) = (X(1) - X(0)) / t(1): u(1) = (X(2) - X(1)) / t(2)

u(2) = (u(1) - u(0)) / (t(2) + t(1))

tx1 = u(0) - u(2) * t(1)

u(0) = (Y(1) - Y(0)) / t(1): u(1) = (Y(2) - Y(1)) / t(2)

u(2) = (u(1) - u(0)) / (t(2) + t(1))

ty1 = u(0) - u(2) * t(1)

u(0) = (X(n) - X(n - 1)) / t(n): u(1) = (X(n - 1) - X(n - 2)) / t(n - 1)

u(2) = (u(0) - u(1)) / (t(n) + t(n - 1))

tx2 = u(0) + u(2) * t(n)

u(0) = (Y(n) - Y(n - 1)) / t(n): u(1) = (Y(n - 1) - Y(n - 2)) / t(n - 1)

u(2) = (u(0) - u(1)) / (t(n) + t(n - 1))

ty2 = u(0) + u(2) * t(n)

Case 1 '夹持条件

a(0) = 1: c(0) = 0: dX(0) = tx1: dY(0) = ty1

a(n) = 1: b(n) = 0: dX(n) = tx2: dY(n) = ty2

Case 2 '自由条件

a(0) = 2: c(0) = 1

dX(0) = 3 * (X(1) - X(0)) / t(1): dY(0) = 3 * (Y(1) - Y(0)) / t(1)

a(n) = 2: b(n) = 1

dX(n) = 3 * (X(n) - X(n - 1)) / t(n): dY(n) = 3 * (Y(n) - Y(n - 1)) / t(n)

Case 3 '循环条件

a(0) = 2: c(0) = 1

dX(0) = 3 * (X(1) - X(0)) / t(1) - (t(1) * (X(2) - X(1)) / t(2) - X(1) + X(0)) / (t(1) + t(2))

dY(0) = 3 * (Y(1) - Y(0)) / t(1) - (t(1) * (Y(2) - Y(1)) / t(2) - Y(1) + Y(0)) / (t(1) + t(2))

a(n) = 2: b(n) = 1

dX(n) = 3 * (X(n) - X(n - 1)) / t(n)

dX(n) = dX(n) + (X(n) - X(n - 1) - t(n) * (X(n - 1) - X(n - 2)) / t(n - 1)) / (t(n) + t(n - 1))

dY(n) = 3 * (Y(n) - Y(n - 1)) / t(n)

dY(n) = dY(n) + (Y(n) - Y(n - 1) - t(n) * (Y(n - 1) - Y(n - 2)) / t(n - 1)) / (t(n) + t(n - 1))

End Select

'计算方程组系数阵和常数阵

For i = 1 To n - 1

a(i) = 2 * (t(i) + t(i + 1)): b(i) = t(i + 1): c(i) = t(i)

dX(i) = 3 * (t(i) * (X(i + 1) - X(i)) / t(i + 1) + t(i + 1) * (X(i) - X(i - 1)) / t(i))

dY(i) = 3 * (t(i) * (Y(i + 1) - Y(i)) / t(i + 1) + t(i + 1) * (Y(i) - Y(i - 1)) / t(i))

Next i

'采用追赶法解方程组

c(0) = c(0) / a(0)

For i = 1 To n - 1

a(i) = a(i) - b(i) * c(i - 1): c(i) = c(i) / a(i)

Next i

a(n) = a(n) - b(n) * c(i - 1)

qx(0) = dX(0) / a(0): qy(0) = dY(0) / a(0)

For i = 1 To n

qx(i) = (dX(i) - b(i) * qx(i - 1)) / a(i)

qy(i) = (dY(i) - b(i) * qy(i - 1)) / a(i)

Next i

px(n) = qx(n): py(n) = qy(n)

For i = n - 1 To 0 Step -1

px(i) = qx(i) - c(i) * px(i + 1)

py(i) = qy(i) - c(i) * py(i + 1)

Next i

'计算曲线上点的坐标

For i = 0 To n - 1

bx3 = (3 * (X(i + 1) - X(i)) / t(i + 1) - 2 * px(i) - px(i + 1)) / t(i + 1)

bx4 = ((2 * (X(i) - X(i + 1)) / t(i + 1) + px(i) + px(i + 1)) / t(i + 1)) / t(i + 1) by3 = (3 * (Y(i + 1) - Y(i)) / t(i + 1) - 2 * py(i) - py(i + 1)) / t(i + 1)

by4 = ((2 * (Y(i) - Y(i + 1)) / t(i + 1) + py(i) + py(i + 1)) / t(i + 1)) / t(i + 1) tt = 0

While (tt <= t(i + 1))

cx = X(i) + (px(i) + (bx3 + bx4 * tt) * tt) * tt

cy = Y(i) + (py(i) + (by3 + by4 * tt) * tt) * tt

u1(num) = cx: v1(num) = cy: num = num + 1: tt = tt + 0.5

Wend

u1(num) = X(i + 1): v1(num) = Y(i + 1): num = num + 1

Next i

End Sub

Private Sub Form_Load()

Dim i As Integer

t = 30

ToInit = False

'Picture1.Scale (0, 0)-(640, 550)

Randomize Timer

https://www.wendangku.net/doc/7d2289562.html,mand2.Caption = "开始初始化"

For i = 0 To t

X(i) = Rnd(1) * 500 + Rnd(1) * 50 + 12

Y(i) = Rnd(1) * 400 + Rnd(1) * 100 + 12

Next i

For i = 0 To t

X(i) = i * 30 + 20

Y(i) = i * 20 + 20

Next i

'Me.Picture1.Picture = LoadPicture("c:\my documents\MenuBack.bmp") Me.Picture1.BackColor = QBColor(0)

End Sub

Private Sub Form_Resize()

On Error Resume Next

Me.Picture1.Height = Me.ScaleHeight - 40

End Sub

Private Sub Form_Unload(Cancel As Integer)

End

End Sub

Private Sub Picture2_MouseDown(Index As Integer, Button As Integer, Shift As Integer, X As Single, Y As Single)

On Error Resume Next

If Button = 1 Then

DownX = X

DownY = Y

Picture2(Index).ZOrder 0

Picture2(Index - 1).BackColor = QBColor(12)

Picture2(Index + 1).BackColor = QBColor(12)

lblX.Caption = "X: " + CStr(CInt(Picture2(Index).Left + 4))

lblY.Caption = "Y: " + CStr(CInt(Picture2(Index).Top + 4))

lblX.Refresh

lblY.Refresh

MovePic Index

Drawposi Index

End If

End Sub

Private Sub Picture2_MouseMove(Index As Integer, Button As Integer, Shift As Integer, X As Single, Y As Single)

If Button = 1 Then

Picture2(Index).Left = Picture2(Index).Left - DownX + X

Picture2(Index).Top = Picture2(Index).Top - DownY + Y

MovePic Index

Command1_Click

Drawposi Index

End If

End Sub

Private Sub Picture2_MouseUp(Index As Integer, Button As Integer, Shift As Integer, X As Single, Y As Single)

On Error Resume Next

If Button = 1 Then

DownX = X

DownY = Y

Picture2(Index - 1).BackColor = QBColor(15)

Picture2(Index + 1).BackColor = QBColor(15)

'MovePic Index

lblX.Caption = "X:"

lblY.Caption = "Y:"

lblX.Refresh

lblY.Refresh

Command1_Click

End If

End Sub

三次样条插值函数

高次插值函数的计算量大,有剧烈振荡,且数值稳定性差；在分段插值中，分段线性插值在分段点上仅连续而不可导，分段三次埃尔米特插值有连续的一阶导数，如此光滑程度常不能满足物理问题的需要，样条函数可以同时解决这两个问题,使插值函数既是低阶分段函数,又是光滑的函数，并且只需在区间端点提供某些导数信息。

5.6.1三次样条函数

定义设在区间［a,b］上取 n+1 个节点

a= =b

函数y=f(x)在各个节点处的函数值为 =f( )(i=0，1，…，n)，若S(x)满足：

(1) S( )= ，i=0，1，…，n ；

(2) 在区间［a,b］上，S(x)具有连续的二阶导数；

(3) 在区间［］(i=0，1，…，n-1)上,S(x) 是x的三次多项式；

称S(x) 是函数y=f(x)在上［a,b］上的三次样条插

S′( +0)=S′( -0)，（5）

S″( +0)=S″( -0)

此时，这样确定的样条函数S(x)称为周期样条函数。

5.6.2三次样条函数的计算

1 用节点处的二阶导数表示的三次样条函数——三弯矩方程

由于S(x)的二阶导数连续，设S(x)在节点的二阶导数为，即

S″( )＝，i=0，1，…，n，

是未知、待定的数。因S(x)是分段三次多项式，则S″(x)是分段一次多项式，且在每个区间［］上，S″(x)可表示为

记＝ - ，则

将上式在区间［］上积分两次，并且由S( )＝，S( ) ＝ 来确定两个积分常数，可得当x∈［］时,

S(x)= (5.6.6)

对上式求导得:

S′(x)= （7）

利用S(x)一阶导数连续的性质，在上式中令x= ，得:

将式（7）中的i换成i-1，得S′(x)在［］上的表达式

S′(x)=

用x= 代入，得：

利用S′( －0)＝S′( ＋0)可得:

两边乘以，得:

i=1,2,…,n-1 (5.6.8)

其中

(9)

这是含有n+1个未知量共有n-1个方程组成的线性方程组，欲确定该方程组的解,尚缺2个方程。因此求三次样条函数还要2个附加条件。

常见的问题有下面两种提法:

①第一类问题附加条件①，即给出边界端点的一阶导数值:S′( )＝，S′( 。

利用前面已推导的公式，当x∈［］时，

S′(x)=

取i＝0，x= ，得:

＝

取i＝n-1，x＝，得：

移项，得：

其中与式(8)联立得一n+1元线性方程组：（10）

其系数矩阵是严格对角占优的三对角矩阵：

可以用追赶法解出(i=0,1,…,n)，将其代入式(6)，即得到三次样条函数的分段表达式。

②第二类问题附加条件②，即给出边界端点的二阶导数值：S″( ) ＝，S″( )＝，代入式(8)得一n-1元线性方程组：

（11）

其系数矩阵为

这是一个三对角矩阵，由于，因而它是严格对角占优的。原方程组是个三对角方程组，可以用追赶法解出，代入式(6)，即得到三次样条函数的分段表达式。

2 用节点处的一阶导数表示的三次样条函数——三转角方程

由于S(x)的一阶导数连续，设S(x)在节点处的一阶导数值为，即

S′( )＝，i＝0，1，…，n，

是未知、待定的数。因S(x)是分段三次多项式，则在每个区间［］上是三次多项式，且满足

S( )＝，S( )＝，S′( )＝，S′( )＝，

故S(x)是［ , ］上的分段三次埃尔米特插值多项式，当x∈［］时

S(x)= （） + ( ) +

(x- ) + (x- ) (i=0，1，…，n-1) (3) 或写为

S(x)= +

将上式在区间［］上求导两次,可得当x∈［］时,

S′(x)= +

和

S″(x) = +

利用S(x)二阶导数连续的性质，在上式中令x= ，得：

S″( +0)=

将上式中的i换成i-1 ，得S″(x)在［］上的表达式

S″(x) == +

用x= 代入，得：

S″( -0)=

利用S″( -0)=S″( +0)可得：

两边乘以，得：

(13)

其中

(14)

这是含有n+1个未知量共有n-1个方程组成的线性方程组，欲确定该方程组的解，尚缺2个方程。因此，求三次样条函数还要2个附加条件。有关情形与节点处的二阶导数表示的三次样条函数类似，叙述如下：

三次样条插值---matlab实现

计算方法实验—三次样条插值 机电学院075094-19 苏建加 20091002764 题目：求压紧三次样条曲线，经过点（-3,2），（-2,0），（1,3），（4,1），而且一阶导 数边界条件S'(-3)=-1;S'（4）=1。 解：首先计算下面的值： 记 1--=j j j x x h ； 1++=j j j j h h h u ；1=+j j u λ ； ?? ????????---+=-++++-j j j j j j j j j j j h y y h y y h h x x x f 1111 111],,[ ；M j =)(''j x s ;],,[611+-=j j j j x x x f d ; h1=-2-(-3)=1;h2=1-(-2)=3;h3=4-1=3; u1=1/4;u2=3/6; d1=6/4*(3/3-(-2)/1)=4.5;d2=6/6*(-2/3-3/3)=-5/3; 由于边界条件S'(-3)=-1;S'（4）=1,得到如下 式子： d0=6/1*(-2/1-(-1))=-6； d3=6/3*(1-(-2)/3)=10/3; 所以得到4个含参数m0~m3 的线性代数方程组为： 2.0000 1.0000 0 0 m0 0.2500 2.0000 0.7500 0 m1 0 0.5000 2.0000 0.5000 m2 0 0 1.0000 2.0000 m3 利用matlab 求解方程得： m = -4.9032 3.8065 -2.5161 2.9247 所以 S1(x)=-0.8172*(-2-x)^3+ 0.6344*(x+3)^3+2.8172*(-2-x)-0.6344*(x+3) x ∈[-3，-2] S2(x)=0.2115*(1-x)^3 -0.1398*(x+2)^3- 1.9032*(1-x)+ 2.2581*(x+2) x ∈[-2，1] S3(x)=-0.1398*(4-x)^3+0.1625(x-1)^3+ 2.2581*(4-x)-1.1290*(x-1) x ∈[1，4] 化简后得：S1(x)=1.4516*x^3 + 10.6128*x^2 + 23.4836*x + 16.1288 x ∈[-3，-2] S2(x)=-0.3513x^3-0.2043x^2+1.8492x+1.7061 x ∈[-2，1] S3(x)=0.3023x^3-2.1651x^2+3.8108x+1.0517 x ∈[1，4] 画图验证：

样条插值函数与应用

样条插值函数及应用

摘要 样条函数具有广泛的应用，是现代函数论的一个十分活跃的分支，是计算方法的主要基础和工具之一，由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。 在实际问题中所遇到许多函数往往很复杂，有些甚至是很难找到解析表达式的。根据函数已有的数据来计算函数在一些新的点处的函数值，就是插值法所需要解决的问题。 插值法是数值逼近的重要方法之一，它是根据给定的自变量值和函数值,求取未知函数的近似值。早在一千多年前,我国科学家就在研究历法时就用到了线性插值和二次插值。而在实际问题中,有许多插值函数的曲线要求具有较高的光滑性,在整个曲线中,曲线不但不能有拐点,而且曲率也不能有突变。因此，对于插值函数必须二次连续可微且不变号 ,这就需要用到三次样条插值。 关键词三次样条函数；插值法

目录 引言 0 第一章三次样条插值 (1) 1.1 样条插值函数简介 (1) 1.2 三次样条函数应用 (2) 第二章AMCM91A 估计水塔水流量 (4) 2.1 理论分析及计算 (5) 2.2运用MATLAB软件计算 (8) 参考文献 (13)

引言 样条函数具有广泛的应用，是现代函数论的一个十分活跃的分支，是计算方法的主要基础和工具之一，由于生产和科学技术向前发展的推动以及电子计算机广泛应用的需要，人们便更多地应用这个工具，也更深刻的认识了它的本质。上世纪四十年代，在研究数据处理的问题中引出了样条函数，例如，在1946年Schoenberg将样条引入数学,即所谓的样条函数，直到五十年代，还多应用于统计数据的处理方面，从六十年代起，在航空、造船、汽车等行业中，开始大量采用样条函数。 在我国，从六十年代末开始，从船体数学放样到飞机外形设计，逐渐出现了一个使用样，逐渐出现了一个使用样条函数的热潮，并推广到数据处理的许多问题中。 在实际生活中有许多计算问题对插值函数的光滑性有较高的要求，例如飞机机翼外形、发动机进、排气口都要求有连续的二阶导数，用三次样条绘制的曲线不仅有很好的光滑度，而且当节点逐渐加密时其函数值整体上能很好地逼近被插函数，相应的导数值也收敛于被插函数的导数值，不会发生“龙格现象”。 现在国内外学者对这方面的研究也越来越重视，根据我们的需要来解决不同的问题，而且函数的形式也在不断地改进，长期以来很多学者致力于样条插值的研究，对三次样条的研究已相当成熟。

对样条函数及其插值问题的一点认识

对样条函数及其插值问题的一点认识 样条函数是计算数学以及计算机辅助设计几何设计的重要工具。1946年，I. J. Schoenberg 著名的关于一元样条函数的奠定性论文“Contribution to the problem of application of equidistant data by analytic functions ”发表，建立了一元样条函数的理论基础。自此以后，关于样条函数的研究工作逐渐深入。随着电子计算机技术的不断进步，样条函数的理论以及应用研究得到迅速的发展和广泛的应用。经过数学工作者的努力，已经形成了较为系统的理论体系。 所谓（多项式）样条函数，乃指具有一定光滑性的分段（分片）多项式。一元n 次且n －1阶连续可微的样条函数具有如下的表示式： 1()()()()N n n j j j s x p x c x x x +==+--∞<<+∞∑[] 011,00,01,,...,,(1),...,(),,...,,n n n n N n N N u un u u u u x x x x x S x x x x ++++ +≥??=??

Matlab中插值函数汇总和使用说明.

告: Matlab中插值函数汇总和使用说明收藏 命令1 interp1 功能一维数据插值(表格查找。该命令对数据点之间计算内插值。它找出一元函数f(x在中间点的数值。其中函数f(x由所给数据决定。x:原始数据点 Y:原始数据点 xi:插值点 Yi:插值点 格式 (1yi = interp1(x,Y,xi 返回插值向量yi,每一元素对应于参量xi,同时由向量x 与Y 的内插值决定。参量x 指定数据Y 的点。 若Y 为一矩阵,则按Y 的每列计算。yi 是阶数为length(xi*size(Y,2的输出矩阵。 (2yi = interp1(Y,xi 假定x=1:N,其中N 为向量Y 的长度,或者为矩阵Y 的行数。 (3yi = interp1(x,Y,xi,method 用指定的算法计算插值: ’nearest’:最近邻点插值,直接完成计算; ’linear’:线性插值(缺省方式,直接完成计算;

’spline’:三次样条函数插值。对于该方法,命令interp1 调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函 数。命令spline 用它们执行三次样条函数插值; ’pchip’:分段三次Hermite 插值。对于该方法,命令interp1 调用函数p chip,用于对向量x 与y 执行分段三次内插值。该方法保留单调性与数据的外形; ’cubic’:与’pchip’操作相同; ’v5cubic’:在MATLAB 5.0 中的三次插值。 对于超出x 范围的xi 的分量,使用方法’nearest’、’linear’、’v5cubic’的插值算法,相应地将返回NaN。对其他的方法,interp1 将对超出的分量执行外插值算法。 (4yi = interp1(x,Y,xi,method,'extrap' 对于超出x 范围的xi 中的分量将执行特殊的外插值法extrap。 (5yi = interp1(x,Y,xi,method,extrapval 确定超出x 范围的xi 中的分量的外插值extrapval,其值通常取NaN 或0。 例1 1.>>x = 0:10; y = x.*sin(x; 2.>>xx = 0:.25:10; yy = interp1(x,y,xx; 3.>>plot(x,y,'kd',xx,yy 复制代码 例2 1.>> year = 1900:10:2010;

三次样条插值课后题集

例1 设)(x f 为定义在[0,3]上的函数，有下列函数值表： 且2.0)('0=x f ，1)('3-=x f ，试求区间[0,3]上满足上述条件的三次样条插值函数)(x s 本算法求解出的三次样条插值函数将写成三弯矩方程的形式： ) ()6()() 6()(6)(6)(211123 13 1j j j j j j j j j j j j j j j j x x h h M y x x h h M y x x h M x x h M x s -- + -- + -+ -= +++++其中，方程中的系数 j j h M 6， j j h M 61+， j j j j h h M y )6(2- ， j j j j h h M y ) 6(211++- 将由Matlab 代码中的变量Coefs_1、Coefs_2、Coefs_3以及Coefs_4的值求出。 以下为Matlab 代码： %============================= % 本段代码解决作业题的例1 %============================= clear all clc % 自变量x 与因变量y ，两个边界条件的取值 IndVar = [0, 1, 2, 3]; DepVar = [0, 0.5, 2, 1.5];

LeftBoun = 0.2; RightBoun = -1; % 区间长度向量，其各元素为自变量各段的长度h = zeros(1, length(IndVar) - 1); for i = 1 : length(IndVar) - 1 h(i) = IndVar(i + 1) - IndVar(i); end % 为向量μ赋值 mu = zeros(1, length(h)); for i = 1 : length(mu) - 1 mu(i) = h(i) / (h(i) + h(i + 1)); end mu(i + 1) = 1; % 为向量λ赋值 lambda = zeros(1, length(h)); lambda(1) = 1; for i = 2 : length(lambda) lambda(i) = h(i) / (h(i - 1) + h(i));

三次样条插值方法的应用

CENTRAL SOUTH UNIVERSITY 数值分析实验报告

三次样条插值方法的应用 一、问题背景 分段低次插值函数往往具有很好的收敛性，计算过程简单，稳定性好，并且易于在在电子计算机上实现，但其光滑性较差，对于像高速飞机的机翼形线船体放样等型值线往往要求具有二阶光滑度，即有二阶连续导数，早期工程师制图时，把富有弹性的细长木条（即所谓的样条）用压铁固定在样点上，在其他地方让他自由弯曲，然后沿木条画下曲线，称为样条曲线。样条曲线实际上是由分段三次曲线并接而成，在连接点即样点上要求二阶导数连续，从数学上加以概括就得到数学样条这一概念。下面我们讨论最常用的三次样条函数及其应用。 二、数学模型 样条函数可以给出光滑的插值曲线（面），因此在数值逼近、常微分方程和偏微分方程的数值解及科学和工程的计算中起着重要的作用。 设区间[]b ,a 上给定有关划分b x x n =<<<= 10x a ，S 为[]b ,a 上满足下面条件的函数。 ● )(b a C S ,2∈； ● S 在每个子区间[]1,+i i x x 上是三次多项式。 则称S 为关于划分的三次样条函数。常用的三次样条函数的边界条件有三种类型： ● Ⅰ型 ()()n n n f x S f x S ''0'',==。 ● Ⅱ型 ()()n n n f x S f x S ''''0'''',==，其特殊情况为()()0''''==n n x S x S 。 ● Ⅲ型 ()() 3,2,1,0,0==j x S x S n j j ，此条件称为周期样条函数。 鉴于Ⅱ型三次样条插值函数在实际应用中的重要地位，在此主要对它进行详细介绍。 三、算法及流程 按照传统的编程方法，可将公式直接转换为MATLAB 可是别的语言即可；另一种是运用矩阵运算，发挥MATLAB 在矩阵运算上的优势。两种方法都可以方便地得到结果。方法二更直观，但计算系数时要特别注意。这里计算的是方法一的程序，采用的是Ⅱ型边界条件，取名为spline2.m 。 Matlab 代码如下： function s=spline2(x0,y0,y21,y2n,x) %s=spline2(x0,y0,y21,y2n,x) %x0,y0 are existed points,x are insert points,y21,y2n are the second

数值分析作业-三次样条插值

数值计算方法作业 实验4.3 三次样条差值函数 实验目的： 掌握三次样条插值函数的三弯矩方法。 实验函数: dt e x f x t ? ∞ -- = 2 221)(π x 0.0 0.1 0.2 0.3 0.4 F(x) 0.5000 0.5398 0.5793 0.6179 0.7554 求f(0.13)和f(0.36)的近似值 实验内容： （1） 编程实现求三次样条插值函数的算法，分别考虑不同的边界条件； （2） 计算各插值节点的弯矩值； （3） 在同一坐标系中绘制函数f(x)，插值多项式，三次样条插值多项式的曲线 比较插值结果。 实验4.5 三次样条差值函数的收敛性 实验目的： 多项式插值不一定是收敛的，即插值的节点多，效果不一定好。对三次样条插值函数如何呢？理论上证明三次样条插值函数的收敛性是比较困难的，通过本实验可以证明这一理论结果。 实验内容： 按照一定的规则分别选择等距或非等距的插值节点，并不断增加插值节点的个数。 实验要求： （1） 随着节点个数的增加，比较被逼近函数和三样条插值函数的误差变化情 况，分析所得结果并与拉格朗日插值多项式比较； （2） 三次样条插值函数的思想最早产生于工业部门。作为工业应用的例子，考 实验名称 实验 4.3三次样条插值函数（P126） 4.5三次样条插值函数的收敛性（P127） 实验时间 姓名 班级 学号 成绩

虑如下例子：某汽车制造商根据三次样条插值函数设计车门曲线，其中一 段数据如下： k x 0 1 2 3 4 5 6 7 8 9 10 k y 0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29 k y ' 0.8 0.2 算法描述： 拉格朗日插值： 错误！未找到引用源。 其中错误！未找到引用源。是拉格朗日基函数，其表达式为：() ∏ ≠=--=n i j j j i j i x x x x x l 0) ()( 牛顿插值： ) )...()(](,...,,[....))(0](,,[)0](,[)()(1102101210100----++--+-+=n n n x x x x x x x x x x f x x x x x x x f x x x x f x f x N 其中????? ?? ?? ?????? --=--= --= -)/(]),...,[],...,[(]...,[..],[],[],,[)()(],[01102110x x x x x f x x x f x x x f x x x x f x x f x x x f x x x f x f x x f n n n n i k j i k j k j i j i j i j i 三样条插值： 所谓三次样条插值多项式Sn(x)是一种分段函数，它在节点Xi(a

关于三次样条插值函数的学习报告(研究生)资料

学习报告—— 三次样条函数插值问题的讨论 班级：数学二班 学号：152111033 姓名：刘楠楠

样条函数： 由一些按照某种光滑条件分段拼接起来的多项式组成的函数；最常用的样条函数为三次样条函数，即由三次多项式组成，满足处处有二阶连续导数。 一、三次样条函数的定义： 对插值区间[,]a b 进行划分，设节点011n n a x x x x b -=<< <<=，若 函数2()[,]s x c a b ∈在每个小区间1[,]i i x x +上是三次多项式，则称其为三次样条函数。如果同时满足()()i i s x f x = (0,1,2)i n =，则称()s x 为()f x 在 [,]a b 上的三次样条函数。 二、三次样条函数的确定： 由定义可设：101212 1(),[,] (),[,]()(),[,] n n n s x x x x s x x x x s x s x x x x -∈??∈?=???∈?其中()k s x 为1[,]k k x x -上的三次 多项式，且满足11(),()k k k k k k s x y s x y --== (1,2,,k n = 由2()[,]s x C a b ∈可得：''''''()(),()(),k k k k s x s x s x s x -+-+== 有''1()(),k k k k s x s x -++= ''''1()(),(1 ,2,,1)k k k k s x s x k n -+ +==-， 已知每个()k s x 均为三次多项式，有四个待定系数，所以共有4n 个待定系数，需要4n 个方程才能求解。前面已经得到22(1)42n n n +-=-个方程，因此要唯一确定三次插值函数，还要附加2个条件，一般上，实际问题通常对样条函数在端点处的状态有要求，即所谓的边界条件。 1、第一类边界条件：给定函数在端点处的一阶导数，即 ''''00(),()n n s x f s x f == 2、第二类边界条件：给定函数在端点处的二阶导数，即

(精选)三次样条插值的MATLAB实现

MATLAB 程序设计期中考查 在许多问题中，通常根据实验、观测或经验得到的函数表或离散点上的信息，去研究分析函数的有关特性。其中插值法是一种最基本的方法，以下给出最基本的插值问题——三次样条插值的基本提法： 对插值区间[]b a ,进行划分：b x x x a n ≤

三次样条插值

三次样条插值 C++数值算法（第二版) 3.3 三次样条插值 给定一个列表显示的函数 yi=y(xi),i=0,1,2,...,N-1。特别注意在xj和xj+1之间的一个特殊的区间。该区间的线性插值公式为(3.3.1)式和(3.3.2)式是拉格朗日插值公式(3.1.1)的特殊情况。 因为它是(分段)线性的，(3.3.1)式在每一区间内的二阶导数为零，在横坐标为xj处的二阶导数不定义或无限。三次样条插值的目的就是要得到一个内插公式，不论在区间内亦或其边界上，其一阶导数平滑，二阶导数连续。 做一个与事实相反的个假设，除yi的列表值之外，我们还有函数二阶导数y"的列表值，即一系列的yi"值，则在每个区间内，可以在(3.3.1)式的右边加上一个三次多项式，其二阶导数从左边的yj"值线性变化到右边的yj+1"值，这么做便得到了所需的连续二阶导数。如果还将三次多项式构造在xj和xj+1处为零，则不会破坏在终点xj和xj+1处与列表函数值yj和yj+1的一致性。 进行一些辅助计算便可知，仅有一种办法才能进行这种构造，即用

注意,(3.3.3)式和(3.3.4)式对自变量x的依赖，是完全通过A和B对x的线性依赖，以及C和D(通过A和B) 对x的三次依赖而实现。 可以很容易地验证，y"事实上是该插值多项式的二阶导数。使用ABCD的定义对x求(3.3.3)式的导数，计算dA/dx dB/dx dC/dx dD/dx,结果为一阶导数因为x=xj是A=1， x=x(i+1)时A=0，而B正相反，则(3.3.6)式表明y"恰为列表函数的二阶导数。而且该二阶导数在两个区间(xj-1, xj)和(xj, xj+1)上是连续的。 现在唯一的问题是，假设yj"是已知的。而实际上并不知道。然而，仍不要求从(3.3.5)式算出的一阶导数在两个区间的边界处是连续的。三次样条的关键思想就在于要求这种连续性，并用它求得等式的二阶导数yi"。 设(3.3.5)式在区间(xj-1, xj)上对x=xj求得的值，等于同一等式在区间(xj,x[j+1])上对x=xj求得的值，便可得到所求方程，是新整理得到(对j=1,......,N-2) 这意味着有N-2个线性方程，但却有N个未知数yi",i=0,......,N-1。因此，存在一个具有两个参数的可能解集。为求得唯一解，需要给出两个进一步的条件，一般取x0和 x[n-1]处的边界条件。最常见的做法有： [1]设y[0"]和y[n-1]"之一或两个都为零，得到所谓的

三次样条插值自然边界条件

例：已知一组数据点，编写一程序求解三次样条插值函数满足 并针对下面一组具体实验数据 0.25 0.3 0.39 0.45 0.53 0.5000 0.5477 0.6245 0.6708 0.7280 求解，其中边界条件为. 1）三次样条插值自然边界条件源程序： function s=spline3(x,y,dy1,dyn) %x为节点，y为节点函数值，dy1，dyn分别为x=0.25,0.53处的二阶导 m=length(x);n=length(y); if m~=n error('x or y输入有误') return end h=zeros(1,n-1); h(n-1)=x(n)-x(n-1); for k=1:n-2 h(k)=x(k+1)-x(k); v(k)=h(k+1)/(h(k+1)+h(k)); u(k)=1-v(k); end g(1)=3*(y(2)-y(1))/h(1)-h(1)/2*dy1; g(n)=3*(y(n)-y(n-1))/h(n-1)+h(n-1)/2*dyn; for i=2:n-1 g(i)=3*(u(i-1)*(y(i+1)-y(i))/h(i)+v(i-1)*(y(i)-y(i-1))/h(i-1)); end for i=2:n-1; A(i,i-1)=v(i-1); A(i,i+1)=u(i-1); end A(n,n-1)=1; A(1,2)=1; A=A+2*eye(n); M=zhuigf(A,g); %调用函数，追赶法求M fprintf('三次样条（三对角）插值的函数表达式\n'); syms X;

for k=1:n-1 fprintf('S%d--%d:\n',k,k+1); s(k)=(h(k)+2*(X-x(k)))./h(k).^3.*(X-x(k+1)).^2.*y(k)... +(h(k)-2*(X-x(k+1)))./h(k).^3.*(X-x(k)).^2.*y(k+1)... +(X-x(k)).*(X-x(k+1)).^2./h(k).^2*M(k)+(X-x(k+1)).*... (X-x(k)).^2./h(k).^2*M(k+1); end s=s.'; s=vpa(s,4); %画三次样条插值函数图像 for i=1:n-1 X=x(i):0.01:x(i+1); st=(h(i)+2*(X-x(i)))./(h(i)^3).*(X-x(i+1)).^2.*y(i)... +(h(i)-2.*(X-x(i+1)))./(h(i)^3).*(X-x(i)).^2.*y(i+1)... +(X-x(i)).*(X-x(i+1)).^2./h(i)^2*M(i)+(X-x(i+1)).*... (X-x(i)).^2./h(i)^2*M(i+1); plot(x,y,'o',X,st); hold on End plot(x,y); grid on %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %调用的函数： %追赶法 function M=zhuigf(A,g) n=length(A); L=eye(n); U=zeros(n); for i=1:n-1 U(i,i+1)=A(i,i+1); end U(1,1)=A(1,1); for i=2:n L(i,i-1)=A(i,i-1)/U(i-1,i-1); U(i,i)=A(i,i)-L(i,i-1)*A(i-1,i); end Y(1)=g(1); for i=2:n Y(i)=g(i)-L(i,i-1)*Y(i-1); end M(n)=Y(n)/U(n,n); for i=n-1:-1:1 M(i)=(Y(i)-A(i,i+1)*M(i+1))/U(i,i);

数值分析实验报告-插值、三次样条(教育教学)

实验报告：牛顿差值多项式&三次样条 问题：在区间[-1,1]上分别取n=10、20用两组等距节点对龙格函数2 1()25f x x 作多项式插值及三次样条插值，对每个n 值，分别画出插值函数及()f x 的图形。 实验目的：通过编程实现牛顿插值方法和三次样条方法，加深对多项式插值的理解。应用所编程序解决实际算例。 实验要求： 1． 认真分析问题，深刻理解相关理论知识并能熟练应用； 2． 编写相关程序并进行实验； 3． 调试程序，得到最终结果； 4． 分析解释实验结果； 5． 按照要求完成实验报告。 实验原理： 详见《数值分析 第5版》第二章相关内容。 实验内容： （1）牛顿插值多项式 1.1 当n=10时： 在Matlab 下编写代码完成计算和画图。结果如下： 代码： clear all clc x1=-1:0.2:1; y1=1./(1+25.*x1.^2); n=length(x1); f=y1(:); for j=2:n for i=n:-1:j f(i)=(f(i)-f(i-1))/(x1(i)-x1(i-j+1)); end end syms F x p ; F(1)=1;p(1)=y1(1); for i=2:n F(i)=F(i-1)*(x-x1(i-1)); p(i)=f(i)*F(i);

end syms P P=sum(p); P10=vpa(expand(P),5); x0=-1:0.001:1; y0=subs(P,x,x0); y2=subs(1/(1+25*x^2),x,x0); plot(x0,y0,x0,y2) grid on xlabel('x') ylabel('y') P10即我们所求的牛顿插值多项式，其结果为：P10（x）=-220.94*x^10+494.91*x^8-9.5065e-14*x^7-381.43*x^6-8.504e-14*x^5+123.36*x^4+2.0202e-1 4*x^3-16.855*x^2-6.6594e-16*x+1.0 并且这里也能得到该牛顿插值多项式的在[-1，1]上的图形，并和原函数进行对比（见Fig.1）。 Fig.1 牛顿插值多项式（n=10）函数和原函数图形 从图形中我们可以明显的观察出插值函数在两端点处发生了剧烈的波动，产生了极大的误差，即龙格现象，当n=20时，这一现象将更加明显。 1.2 当n=20时： 对n=10的代码进行修改就可以得到n=20时的代码。将“x1=-1:0.2:1;”改为“x1=-1:0.1:1;”即可。运行程序，我们得到n=20时的牛顿插值多项式，结果为：P20(x)= 260188.0*x^20 - 1.0121e6*x^18 + 2.6193e-12*x^17 + 1.6392e6*x^16 + 2.248e-11*x^15 - 1.4429e6*x^14 - 4.6331e-11*x^13 + 757299.0*x^12 + 1.7687e-11*x^11 - 245255.0*x^10 + 2.1019e-11*x^9 + 49318.0*x^8 + 3.5903e-12*x^7 - 6119.2*x^6 - 1.5935e-12*x^5 + 470.85*x^4 + 1.3597e-14*x^3 - 24.143*x^2 - 1.738e-14*x + 1.0 同样的，这里得到了该牛顿插值多项式的在[-1，1]上的图形，并和原函数进行对比（见Fig.2）。

Matlab中插值函数汇总和使用说明

告： Matlab中插值函数汇总和使用说明收藏 命令1 interp1 功能一维数据插值（表格查找）。该命令对数据点之间计算内插值。它找出一元函数f(x)在中间点的数值。其中函数f(x)由所给数据决定。x：原始数据点 Y：原始数据点 xi：插值点 Yi：插值点 格式 (1)yi = interp1(x,Y,xi) 返回插值向量yi，每一元素对应于参量xi，同时由向量x 与Y 的内插值决定。参量x 指定数据Y 的点。 若Y 为一矩阵，则按Y 的每列计算。yi 是阶数为length(xi)*size(Y,2)的输出矩阵。 (2)yi = interp1(Y,xi) 假定x=1:N，其中N 为向量Y 的长度，或者为矩阵Y 的行数。 (3)yi = interp1(x,Y,xi,method) 用指定的算法计算插值： ’nearest’：最近邻点插值，直接完成计算； ’linear’：线性插值（缺省方式），直接完成计算； ’spline’：三次样条函数插值。对于该方法，命令interp1 调用函数spline、ppval、mkpp、umkpp。这些命令生成一系列用于分段多项式操作的函

数。命令spline 用它们执行三次样条函数插值； ’pchip’：分段三次Hermite 插值。对于该方法，命令interp1 调用函数p chip，用于对向量x 与y 执行分段三次内插值。该方法保留单调性与数据的外形； ’cubic’：与’pchip’操作相同； ’v5cubic’：在MATLAB 5.0 中的三次插值。 对于超出x 范围的xi 的分量，使用方法’nearest’、’linear’、’v5cubic’的插值算法，相应地将返回NaN。对其他的方法，interp1 将对超出的分量执行外插值算法。 (4)yi = interp1(x,Y,xi,method,'extrap') 对于超出x 范围的xi 中的分量将执行特殊的外插值法extrap。 (5)yi = interp1(x,Y,xi,method,extrapval) 确定超出x 范围的xi 中的分量的外插值extrapval，其值通常取NaN 或0。 例1 1.>>x = 0:10; y = x.*sin(x); 2.>>xx = 0:.25:10; yy = interp1(x,y,xx); 3.>>plot(x,y,'kd',xx,yy) 复制代码 例2 1.>> year = 1900:10:2010; 2.>> product = [75.995 91.972 105.711 12 3.203 131.669 150.697 179.323 203.212 226.505

三次样条插值函数matlab程序绝不坑爹

x0=[0 0.9211 1.8431 2.9497 3.8714 4.9781 5.9 7.0064 7.9286 8.9678 10.9542 12.0328 12.9544 13.8758 14.9822 15.9039 16.8261 17.9317 19.0375 19.9594 20.8392 22.9581 23.88 24.9869 25.9083]; >> >> y0=[14405 11180 10063 11012 8797 9992 8124 10160 8488 11018 19469 20196 18941 15903 18055 15646 13741 14962 16653 14496 14648 15225 15264 13708 9633]; >> x=0:0.1:25.9; >> y1=interp1(x0,y0,x,'spline'); >> pp1=csape(x0,y0); %样条插值工具箱函数 y2=ppval(pp1,x); %计算x对应的y值 pp2=csape(x0,y0,'second'); y3=ppval(pp2,x); xydata=[x',y1',y2',y3'] subplot(1,2,1) plot(x0,y0,'+',x,y1) title('Spline1') subplot(1,2,2) plot(x0,y0,'+',x,y2) title('Spline2') dx=diff(x); dy=diff(y2); dy_dx=dy./dx; dy_dx0=dy_dx(1) ytemp=y2(13<=x&x<=15); ymin=min(ytemp); xmin=x(y2==ymin); xymin_1315=[xmin,ymin]

实验四三次样条插值Word版

实验四三次样条插值的应用 一、问题描述 The upper portion of this noble beast is to be approximated using clamped cubic spline interpolants. The curve is drawn on a grid from which the table is constructed. Use Algorithm 3.5 to construct the three clamped cubic splines. 二、模型建立 三次样条插值 给定一个列表显示的函数 yi=y(xi),i=0,1,2,...,N-1。特别注意在xj和xj+1之间的一个特殊的区间。该区间的线性插值公式为：

(3.3.1)式和(3.3.2)式是拉格朗日插值公式(3.1.1)的特殊情况。 因为它是(分段)线性的，(3.3.1)式在每一区间内的二阶导数为零，在横坐标为xj处的二阶导数不定义或无限。三次样条插值的目的就是要得到一个内插公式，不论在区间内亦或其边界上，其一阶导数平滑，二阶导数连续。 做一个与事实相反的个假设，除yi的列表值之外，我们还有函数二阶导数y"的列表值，即一系列的yi"值，则在每个区间内，可以在(3.3.1)式的右边加上一个三次多项式，其二阶导数从左边的yj"值线性变化到右边的yj+1"值，这么做便得到了所需的连续二阶导数。如果还将三次多项式构造在xj和xj+1处为零，则不会破坏在终点xj和xj+1处与列表函数值yj和yj+1的一致性。 进行一些辅助计算便可知，仅有一种办法才能进行这种构造，即用 注意,(3.3.3)式和(3.3.4)式对自变量x的依赖，是完全通过A和B对x的线性依赖，以及C和D(通过A和B)对x的三次依赖而实现。可以很容易地验证，y"事实上是该插值多项式的二阶导数。使用ABCD的定义对x求(3.3.3)式的导数，计算dA/dx dB/dx dC/dx dD/dx,结果为一阶导数

三次样条插值函数

一．介绍

二．程序框图 开始 输入未知数X及 (xi,yi),i=0,1,…,n 计算步长H[i] 计算λ、μ、d 根据边界条件，求 解相应的方程得到 M1,…, Mn 将M代入原方程， 得到分段函数 结束

三．源码 syms h n=9;%插入节点数，可以根据题目更改 h=2/(n+1); u=0.5; v=0.5; f=inline('1/(1+25*x.^2)');%输入函数，这个也可以根据题目更改g=inline('3/h*((c-b)/h-(b-a)/h)','a','b','c','h'); for i=1:n+2 x(1)=-1; x(i+1)=x(i)+2/(n+1); y(i)=f(x(i)); end for i=1:n d(i)=g(y(i),y(i+1),y(i+2),h); end A=zeros(n,n); for i=1:n A(i,i)=2; end for i=1:n-1 A(i,i+1)=u; A(i+1,i)=v; end B=zeros(n,1); for i=1:n B(i,1)=d(i) end C=inv(A)*B for i=1:n M(i)=C(i,1); end x=(-1:h/50:1); k=1./(1+25*x.^2); cs=spline(x,k); plot(x,k,'r.'); hold on; ezplot('1/(1+25*x^2)',[-1 1]); title('三次样条插值曲线和f(x)曲线') 四．结果

系数矩阵 弯矩M

分段函数不同次幂x对应的系数 三次样条插值函数与原函数

三次样条插值多项式matlab

三次样条插值多项式 ——计算物理实验作业四 陈万物理学2013级 主程序： clear,clc; format rat x = [1,4,9,16,25,36,49,64]; y = [1,2,3,4,5,6,7,8]; f1 = ; fn = 1/16; [a,b,c,d,M,S] = spline(x,y,f1,fn); 子程序1： function [a,b,c,d,M,S]=spline(x,y,f1,fn) % 三次样条插值函数 % x是插值节点的横坐标 % y是插值节点的纵坐标 % u是插值点的横坐标 % f1是左端点的一阶导数 % fn是右端点的一阶导数 % a是三对角矩阵对角线下边一行 % b是三对角矩阵对角线 % c是三对角矩阵对角线上边一行 % S是插值点的纵坐标

n = length(x); h = zeros(1,n-1); deltay = zeros(1,n); miu = zeros(1,n-1); lamda = zeros(1,n-1); d = zeros(1,n-1); for j = 1:n-1 h(j) = x(j+1)-x(j); deltay(j) = y(j+1)-y(j); end % 得到h矩阵 for j = 2:n-1 sumh = h(j-1) + h(j); miu(j) = h(j-1) / sumh; lamda(j) = h(j) / sumh; d(j) = 6*( deltay(j)/h(j)-(deltay(j-1)/h(j-1)))/sumh; end % 根据第一类边界条件，作如下规定 lamda(1) = 1; d(1) = 6*(deltay(1)/h(1)-f1)/h(1); miu(1) = 1; d(n) = 6*(fn-deltay(n-1)/h(n-1))/h(n-1);