基于李庆扬所著《数值分析》。ps由于比较忙，只是编写出来，没有时间找bug和优化。

一、拉格朗日插值

lagrange

#依次为观测值的横坐标，纵坐标，自变量(n2维)取值

n1

result

for(i in 1:n2){

constant1[i]

for(j in

1:n1){

constant2[j]

temp

temp

result[i]

}

result[i]

}

return(result)

}

x

y

inv

lagrange(x,y,inv)

二、牛顿均差插值

niux

#定义函数

k

#定义循环次数为k

A

#定义均差表A

s

#定义最终结果为s

r

#式子中一个因子r

w

#式子中一个因子w

for(j in 1:k+1){

#均差表一共有k+1列

for(i in 1:k){

#均差表一共有k行

A[i,1]

#均差表第一列，第二列分别为x和f(x)

if(j<=2|j>=i+2){A[i,j]

#均差表左上角部分

else

{A[i,j]

#均差表内元素

}

}

for(i in 2:k) {w[i]

#牛顿插值公式的一个因子w

for(i in 1:k) {r[i]

#牛顿插值公式另一个因子r

s

#最终结果s

list(均差表=A,结果=s)

#输出均差表，结果

}

f

#函数值f

x

b

niux(x,f,b)

三、牛顿差分插值

#牛顿差分形式的插值

NewtonDiff

n1

result

for (i in 1:n2){

if (direct==1){

#向前差分

t

difftable

difftable[,1]

for(col in 2:n1){

for(row in

1:n1){

if

(row+col>(n1+1)) break

difftable[row,col]

}

}

ditui

for(j in 2:n1){

ditui[j]

}

result[i]

}

else

{t

difftable

difftable[,1]

for(col in 2:n1){

for(row in

1:n1){

if (row

difftable[row,col]

}

}

ditui

for(j in 2:n1){

ditui[j]

}

result[i]

}

}

return(result)

}

x

y

inv

direct

NewtonDiff(x,y,inv,direct)

四、分段线性插值

pw_li

ld1

result

juzhen

x

for(i in 1:ld2){

if

(inv[i]>x[ld1-1]){

result[i]

}

else if (inv[i]

result[i]

}

else {

for(j in

2:(ld1-2)){

if

(inv[i]

result[i]

}

}

}

}

return(result)

}

x

y

inv

pw_li(x,y,inv)

五、分段埃尔米特插值

#分段三次埃尔米特插值

hermite

#依次为自变量(增序)，因变量，插值点导数，待计算自变量

ld1

result

for(i in 1:ld2){

if

(inv[i]>x[ld1-1]){

result[i]

((inv[i]-x[ld1-1])/(x[ld1]-x[ld1-1]))^2*y[ld1]*(1+2*(inv[i]-x[ld1])/(x[ld1-1]-x[ld1]))+

((inv[i]-x[ld1])/(x[ld1-1]-x[ld1]))^2*(inv[i]-x[ld1-1])*derv[ld1-1]+

((inv[i]-x[ld1-1])/(x[ld1]-x[ld1-1]))^2*(inv[i]-x[ld1])*derv[ld1]

}

else if (inv[i]

result[i]

((inv[i]-x[1])/(x[2]-x[1]))^2*y[2]*(1+2*(inv[i]-x[2])/(x[1]-x[2]))+

((inv[i]-x[2])/(x[1]-x[2]))^2*(inv[i]-x[1])*derv[1]+

((inv[i]-x[1])/(x[2]-x[1]))^2*(inv[i]-x[2])*derv[2]

}

else {

for(j in

2:(ld1-2)){

if

(inv[i]

result[i]

((inv[i]-x[j])/(x[j+1]-x[j]))^2*y[j+1]*(1+2*(inv[i]-x[j+1])/(x[j]-x[j+1]))+

((inv[i]-x[j+1])/(x[j]-x[j+1]))^2*(inv[i]-x[j])*derv[j]+

((inv[i]-x[j])/(x[j+1]-x[j]))^2*(inv[i]-x[j+1])*derv[j+1]

break

}

}

}

}

return(result)

}

x

y

derv

inv

hermite(x,y,derv,inv)

六、三次样条插值

#三次样条插值

cubicspline

#依次为自变量(增序)，因变量，边界条件(1,2,3),边界取值(二行三列矩阵)和研究自变量

#关于边界条件，2行从上往下为首末端点，三列从左至右为函数值，一阶导数值，二阶导数值

n1

if (cond==1){

d

d[1]

d[n1]

for (i in 2:(n1-1)){

d[i]

((y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]))/(x[i+1]-x[i-1])

)

}

#d-vector is finished

juzhen1=juzhen2=juzhen3

for(i in 1:(n1-2)){

juzhen1[i,i]

juzhen2[i,i+1]

juzhen3[i,i+2]

}

juzhen

juzhen

m

list(m=m)

}

if (cond==2){

d

d[1]

d[n1]

for (i in 2:(n1-1)){

d[i]

((y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]))/(x[i+1]-x[i-1])

)

}

#d-vector is finished

juzhen1=juzhen2=juzhen3

for(i in 1:(n1-2)){

juzhen1[i,i]

juzhen2[i,i+1]

juzhen3[i,i+2]

}

juzhen

juzhen

m

list(m=m)

}

if (cond==3){

d

for(i in 1:(n1-1)){

d[i]

(y[2]-y[1])/(x[2]-x[1])-(y[i+1]-y[i])/(x[i+1]-x[i])

)/(x[2]-x[1]+x[i+1]-x[i])

}

#d-vector is finished

juzhen1=juzhen2=juzhen3

for (i in 1:(n1-3)){

juzhen1[i,i]

juzhen2[i,i+1]

juzhen3[i,i+2]

}

juzhen

juzhen

c(2,0.5,rep(0,n1-4),0.5),juzhen,c((x[2]-x[1])/(x[n1]-x[n1-1]+x[2]-x[1]),rep(0,n1-4),1-(x[2]-x[1])/(x[n1]-x[n1-1]+x[2]-x[1]),2)

)

m

list(m=m)

}

result

for(i in 1:n2){

if (inv[i]

result[i]

m[2]*(inv[i]-x[1])^3/(6*(x[2]-x[1]))+

(y[1]-m[1]*(x[2]-x[1])^2/6)*(x[2]-inv[i])/(x[2]-x[1])+

(y[2]-m[2]*(x[2]-x[1])^2/6)*(inv[i]-x[1])/(x[2]-x[1])

}

else if

(inv[i]>x[n1-1]){

result[i]

m[n1]*(inv[i]-x[n1-1])^3/(6*(x[n1]-x[n1-1]))+

(y[n1-1]-m[n1-1]*(x[n1]-x[n1-1])^2/6)*(x[n1]-inv[i])/(x[n1]-x[n1-1])+

(y[n1]-m[n1]*(x[n1]-x[n1-1])^2/6)*(inv[i]-x[n1-1])/(x[n1]-x[n1-1])

}

else {

for(j in

3:(n1-2)){

if

(inv[i]

result[i]

m[j]*(inv[i]-x[j-1])^3/(6*(x[j]-x[j-1]))+

(y[j-1]-m[j-1]*(x[j]-x[j-1])^2/6)*(x[j]-inv[i])/(x[j]-x[j-1])+

(y[j]-m[j]*(x[j]-x[j-1])^2/6)*(inv[i]-x[j-1])/(x[j]-x[j-1])

}

}

}

}

return(result)

}

x

y

cond

condv

inv

cubicspline(x,y,cond,condv,inv)

如若有兴趣或者提建议，欢迎关注微信公众号：R语言的那些事