数値計算 赤坂隆先生 コロナ社

#newtonの定差内挿公式

#関数定義

f=function(x){
  z<-sin(x)+x*exp(-x^2)+1
  return(z)
}

#差分する終点座標を決定
f_points=15;start_points=10

#分割の細かさを決める
sep_f=0.1

#差分前の元データを作成
f_data=data.frame(x=seq(start_points,f_points,sep_f),f=0)

for(j in 1:nrow(f_data)){

  f_data$f[j]=f(f_data$x[j])

}

#差文商

#元データより細かいメッシュを入れる。
sep=0.01

#結果の実際の値と予測値を入れる箱を作る
func_data=data.frame(x=seq(start_points,f_points,sep),f=0)

#差文商の計算(ここから)

difference=X=array(0,dim=c(nrow(f_data),nrow(f_data)))

difference[1,]=f_data$f;X[1,]=f_data$x

for(j in 1:(nrow(f_data)-1)){

  k=j+1

  vector=difference[j,]

  vector=vector[1:(length(vector)-(j-1))]

  zeros=rep(0,j)

  d_f=c()

  for(l in 1:(length(vector)-1)){

   d_f=c(d_f,(vector[l+1]-vector[l])/sep_f)

   } 

  d_f=c(sep_f*d_f,zeros)

  difference[k,]=d_f

  }

#(ここまで)


#差文商の結果をもとに計算をする(ここから);

  f_0=difference[,1]


  for(j in 1:nrow(func_data)){

    f_value=c();x=func_data$x[j]

    for(k in 1:(length(f_0)-1)){

    f_value=c(f_value,prod(x-f_data$x[1:k])*f_0[k+1]*(1/gamma(k+1))*(1/sep_f)^(k)) 

    }

    f_value=sum(f_value,f_0[1])

    func_data$f[j]=f_value 
  }

func_data=func_data %>% mutate(ff=0)

colnames(func_data)=c("x","predict","samples")

for(j in 1:nrow(func_data)){

  func_data$samples[j]=f(func_data$x[j])

}

#(ここまで)

#図のプロット

plot(func_data$x,func_data$samples,col=2,type="o",ylim=c(min(func_data$predict,func_data$samples),max(func_data$predict,func_data$samples)),main="red;samples , green;predict",xlab="x",ylab="value")
par(new=T)
plot(func_data$x,func_data$predict,col=3,type="o",ylim=c(min(func_data$predict,func_data$samples),max(func_data$predict,func_data$samples)),main="red;samples , green;predict",xlab="x",ylab="value")

#重合公式(一変数関数の積分:台形公式)

#積分関数の定義：
f=function(x){

z<-exp(x)

return(z)

}

#１区間単位の大きさ
h=0.01

#積分区間（終点）
end_point=0.4

#積分区間（始点）
start_point=0

#f0,f1,f2,...を作るデータフレーム
f_data=data.frame(x=seq(start_point,end_point,h),f=0)

for(j in 1:nrow(f_data)){

  f_data$f[j]=f(f_data$x[j])

}

#重合公式の計算
integrated_f=c()

for(j in 1:nrow(f_data)){

  if(j==1|j==nrow(f_data)){

  integrated_f=c(integrated_f,f_data$f[j])

  }else{

  integrated_f=c(integrated_f,2*f_data$f[j])  

  }

}

integrated_f=(h/2)*integrated_f

#補正項がない場合の積分値：
print(paste0("区間",start_point,"から",end_point,"までの積分値は",sum(integrated_f)))

#Rのライブラリ使用
print(paste0("実際の値は",integrate(f,start_point,end_point)$value))

#シンプソンの第一公式による積分値の計算(一変数関数の積分)

#関数の定義：
func=function(x){

y<-exp(x)

return(y)

}

#a:積分区間の始点、a+H:積分区間の終点、n:積分区間の分割数
a=0;H=0.4;n=100000

#重合公式の計算
Y=(H/3)*(1/(n-1))*(func(a)+func(a+H));

for(j in 1:(n-2)){

  #奇数番目は4、偶数番目は2になる係数
  coefficient=ifelse(j-floor(j/2)*2==0,4,2)

  Y=Y+((H*coefficient)/(3*(n-1)))*func(a+(j*H)/(n-1))

}

#simpson積分値
print(Y)

#ライブラリによる計算
integrate(func,a,a+H)

#多重積分

#２重積分(台形公式)

x1=0;x2=1

y1=1;y2=2

h=0.25;k=0.5

f=function(x,y){

  z<-x^2+x*y

  return(z)

}

Y=seq(y1,y2,k)

X=seq(x1,x2,h)

lattice1=array(4,dim=c(length(X),length(Y)))

lattice2=array(0,dim=c(length(X),length(Y)))

for(j in 1:length(X)){
  for(i in 1:length(Y)){

  lattice2[j,i]=f(X[j],Y[i]) 

  }
}

lattice1[,1]=2;lattice1[1,]=2;

lattice1[,ncol(lattice1)]=2;lattice1[nrow(lattice1),]=2

lattice1[1,1]=1;lattice1[nrow(lattice1),ncol(lattice1)]=1;

lattice1[1,ncol(lattice1)]=1;lattice1[nrow(lattice1),1]=1


print(paste0("重積分の値は",sum((1/4)*h*k*lattice1*lattice2)))

#rm(list=ls())

#多重積分の計算　sqrt(x)<=y<=1, 0<=x<=1

#f(x)=sin(pi*x^3/2)

#関数の定義：
func=function(x){

y<-sin(pi*x^3/2)

return(y)

}


simpson=function(x){

#a:積分区間の始点、a+H:積分区間の終点、n:積分区間の分割数
n=10000;a=sqrt(x);H=1-sqrt(x)

#重合公式の計算
Y=(H/3)*(1/(n-1))*(func(a)+func(a+H));

for(j in 1:(n-2)){

  #奇数番目は4、偶数番目は2になる係数
  coefficient=ifelse(j-floor(j/2)*2==0,4,2)

  Y=Y+((H*coefficient)/(3*(n-1)))*func(a+(j*H)/(n-1))

}

#simpson積分値
return(Y)

}


#a:積分区間の始点、a+H:積分区間の終点、n:積分区間の分割数
n=10000;a=0;H=1

#重合公式の計算
Y=(H/3)*(1/(n-1))*(simpson(a)+simpson(a+H));

for(j in 1:(n-2)){

  #奇数番目は4、偶数番目は2になる係数
  coefficient=ifelse(j-floor(j/2)*2==0,4,2)

  Y=Y+((H*coefficient)/(3*(n-1)))*simpson(a+(j*H)/(n-1))

}

#多重積分

#２重積分(シンプソンの第一公式)

x1=0;x2=1

y1=1;y2=2

h=0.01;k=0.01

f=function(x,y){

  z<-1/(x+y)

  return(z)

}

Y=seq(y1,y2,k)

X=seq(x1,x2,h)

lattice1=array(4,dim=c(length(X),length(Y)))

lattice2=array(0,dim=c(length(X),length(Y)))


coefficient1=c(1);coefficient2=c(4);coefficient3=c(2)  

for(j in 2:(length(X)-1)){

coefficient1=c(coefficient1,ifelse(j/2==floor(j/2),4,2))  

coefficient2=c(coefficient2,ifelse(j/2==floor(j/2),16,8))  

coefficient3=c(coefficient3,ifelse(j/2==floor(j/2),8,4))  

}

coefficient1=c(coefficient1,1);coefficient2=c(coefficient2,4);coefficient3=c(coefficient3,2)

lattice1=t(lattice1)

for(j in 1:nrow(lattice1)){

if(j==1|j==nrow(lattice1)){

lattice1[j,]=coefficient1  

}else{

if(j/2==floor(j/2)){

lattice1[j,]=coefficient2  

}else{

lattice1[j,]=coefficient3  

}  

}  

}



for(j in 1:length(X)){
  for(i in 1:length(Y)){

  lattice2[j,i]=f(X[j],Y[i]) 

  }
}


lattice1=t(lattice1)

print(paste0("重積分の値は",sum((1/9)*h*k*lattice1*lattice2)))

#多重積分（モンテカルロ法）

f=function(x){

  z<-exp(-(x[1]*x[3]+sqrt(x[1]*x[2]*x[4])+sin(x[2]*x[3])+cos(x[4]*x[5])+x[5]^2))

  return(z)

}

N=10000

h=0.01

#積分領域
x11=0;x12=1
x21=0;x22=1
x31=0;x32=1
x41=0;x42=1
x51=0;x52=1

f_value=c();


for(j in 1:N){

x=array(0,dim=c(1,5))

x[1]=sample(seq(x11,x12,h),1)

x[2]=sample(seq(x21,x22,h),1)

x[3]=sample(seq(x31,x32,h),1)

x[4]=sample(seq(x41,x42,h),1)

x[5]=sample(seq(x51,x52,h),1)

f_value=c(f_value,f(x))


}



alpha=0.95



print(paste0(alpha*100,"%区間は(",(1/N)*sum(f_value)+qnorm((1-alpha)/2)*(1/sqrt(N))*(sqrt(sum((f_value-(1/N)*sum(f_value))^2)/(N-1)))," , ",(1/N)*sum(f_value)+qnorm(alpha+(1-alpha)/2)*(1/sqrt(N))*(sqrt(sum((f_value-(1/N)*sum(f_value))^2)/(N-1))),")"))

print(paste0("推定誤差は",qnorm(alpha+(1-alpha)/2)*(1/sqrt(N))*(sqrt(sum((f_value-(1/N)*sum(f_value))^2)/(N-1)))))

print(paste0("推定値は",(1/N)*sum(f_value)))

#フーリエ近似(データに対する解析関数近似)

library(dplyr)

#サンプルデータに対しフーリエ級数で近似する
data=data.frame(x=seq(-20,20,1),f=0)

f=function(x){

  z<-20*sin(8*x)+9*cos(9*x)+2*cos(10*x)

  return(z)

}

for(j in 1:nrow(data)){

  data$f[j]=f(data$x[j])

}

N=length(data$x[sign(data$x)==1])

A0=(1/(2*N))*sum(data$f)

A=c()

B=c()

for(j in 1:N){

  A=c(A,(1/N)*sum(data$f*cos(j*data$x)))

  B=c(B,(1/N)*sum(data$f*sin(j*data$x)))

}

fourier_trans=function(x){

 z<-A0

 for(j in 1:N){

   z<-z+A[j]*cos(j*x)+B[j]*sin(j*x)
 }

 return(z) 

}

data=data %>% dplyr::mutate(predict=0)

for(j in 1:nrow(data)){

  data$predict[j]=fourier_trans(data$x[j])

}


plot(data$x,data$f,type="l",col=2,xlab="x",ylab="values",main="red:predict, green:samples",ylim=c(min(data$f,data$predict),max(data$f,data$predict)))
par(new=T)
plot(data$x,data$predict,type="l",col=3,xlab="x",ylab="values",main="red:predict, green:samples",ylim=c(min(data$f,data$predict),max(data$f,data$predict)))

print(paste0("予測データと実際のデータとの相関は",cor(data$f,data$predict)))

#F_test(A=B=0)

f_value_data=data.frame(A=A,B=B,f=0,test=0)

alpha=0.1

for(j in 1:N){

A_sub=A;A_sub[j]=0

B_sub=B;B_sub[j]=0

f_value_data$f[j]=(2*N-2)*sum(A[j]^2+B[j]^2)/(2*sum(A_sub^2+B_sub))

f_value_data$test[j]=qf(1-alpha,2,2*N-2) 

}

f_value_data=f_value_data %>% mutate(A_sub=ifelse(f>test,A,0),B_sub=ifelse(f>test,B,0))

A=f_value_data$A_sub

B=f_value_data$B_sub

fourier_trans=function(x){

 z<-A0

 for(j in 1:N){

   z<-z+A[j]*cos(j*x)+B[j]*sin(j*x)
 }

 return(z) 

}

data=data %>% dplyr::mutate(predict=0)

for(j in 1:nrow(data)){

  data$predict[j]=fourier_trans(data$x[j])

}


plot(data$x,data$f,type="l",col=2,xlab="x",ylab="values",main="red:predict, green:samples",ylim=c(min(data$f,data$predict),max(data$f,data$predict)))
par(new=T)
plot(data$x,data$predict,type="l",col=3,xlab="x",ylab="values",main="red:predict, green:samples",ylim=c(min(data$f,data$predict),max(data$f,data$predict)))

#チェビシェフ近似(データに対する解析関数近似)


library(dplyr)

#サンプルデータに対しチェビシェフ多項式で近似する

sep=0.001

#多項式の次数

r=50

data=data.frame(x=seq(-1,1,sep),f=0)


f=function(x){

  z<-x*sin(x)

  return(z)

}

for(j in 1:nrow(data)){

  data$f[j]=f(data$x[j])

}

c0=(1/(nrow(data)+1))*sum(data$f)

c=c()

for(j in 1:r){


T_f=function(x){

  T1=x;T2=2*x^2-1;

  T_vec=c(T1,T2)

  if(j==1){

    T_val=T1

  }else{

    if(j==2){

    T_val=T2

    }else{



  for(k in 1:(j-2)){

    T_val=2*x*T_vec[2]-T_vec[1]

    T_vec=c(T_vec[2],T_val)

  }



}

}

  y<-T_val

  return(y)

}

C_k=c()

for(l in 1:nrow(data)){

  C_k=c(C_k,data$f[l]*T_f(data$x[l]))

}

c=c(c,(2/(nrow(data)+1))*sum(C_k))

}

data=data %>% dplyr::mutate(predict=0)

t_f=array(0,dim=c(nrow(data),r+2))

t_f[,1]=data$x;t_f[,2]=rep(c0,nrow(data))

for(j in 1:r){


T_f=function(x){

  T1=x;T2=2*x^2-1;

  T_vec=c(T1,T2)

  if(j==1){

    T_val=T1

  }else{

    if(j==2){

    T_val=T2

    }else{



  for(k in 1:(j-2)){

    T_val=2*x*T_vec[2]-T_vec[1]

    T_vec=c(T_vec[2],T_val)

  }



}

}

  y<-T_val

  return(y)

}

for(l in 1:nrow(t_f)){

t_f[l,(1+1+j)]=c[j]*T_f(t_f[l,1])

}

}

data=data %>% dplyr::mutate(predict=apply(t_f[,2:(r+2)],1,sum))

plot(data$x,data$f,type="l",col=2,ylim=c(min(data$f,data$predict),max(data$f,data$predict)),xlab="x",ylab="value",main=paste0("Dim=",r,"; red:samples, green:predict"))
par(new=T)
plot(data$x,data$predict,type="l",col=3,ylim=c(min(data$f,data$predict),max(data$f,data$predict)),xlab="x",ylab="value",main=paste0("Dim=",r,"; red:samples, green:predict"))

T_f=function(x){

  T1=x;T2=2*x^2-1;

  T_vec=c(T1,T2)

  if(j==1){

    T_val=T1

  }else{

    if(j==2){

    T_val=T2

    }else{



  for(k in 1:(r+1-2)){

    T_val=2*x*T_vec[2]-T_vec[1]

    T_vec=c(T_vec[2],T_val)

  }



}

}

  y<-T_val

  return(y)

}


#許容誤差限界
error_maximum=abs(T_f(data$x[data$predict==max(data$predict)])[1]*exp(data$x[data$x==max(data$x)]))/(gamma(r+2)*2^(r))

print(paste0("許容誤差限界は",error_maximum))

#ガウス・ザイデル法 p.248

A=matrix(c(10,2,1,1,20,3,2,1,25),ncol=3)

C=c(61.7,200.5,36.2)

r=C/diag(A)

b=A/c(diag(A))

diag(b)=0

mat=cbind(r,-b)

ite=100

x=rep(1,3)

for(l in 1:ite){

z=c(1,x)

x=mat%*%z

print(x)

}

#数値計算　共役傾斜法　p.264

A=matrix(c(5,4,1,4,6,2,1,2,7),ncol=3)

x=c(0,0,0)

c=c(13,17,18)

ite=10

for(l in 1:ite){


if(l==1){

p=c-A%*%x

Ap=A%*%p

a=sum(c^2)

r=c

d=t(p)%*%A%*%p

alpha=as.numeric(a/d)

}else{

x=x+alpha*p

r=r-alpha*A%*%p

a_past=a

a=sum(r^2)

b=a/a_past

p=r+b*p

d=t(p)%*%A%*%p

alpha=as.numeric(a/d)

}

print(x)

}