Go to Qiita Advent Calendar 2024 Top
LoginSignup
1
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@kozakai-ryouta(亮太 小酒井)

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

Last updated at Posted at 2021-04-05
#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)

}
1
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?