r语言安装ipsolve_数值分析的R语言实现（插值部分）
基于李慶揚所著《數值分析》。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語言的那些事
總結
以上是生活随笔為你收集整理的r语言安装ipsolve_数值分析的R语言实现（插值部分）的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： filezilla 设置filezill
- 下一篇： 浅谈python_浅谈python-Dj