流れ解析のための有限要素法

# p.26 修正中です

N=10

x=c(1:N)

u=sample(1:100,N,replace=T)

a_b=array(0,dim=c(2,N-1))

lam=1

integral=c()

for(i in 1:(N-1)){
 
u_vec=u[i:(i+1)];x_vec=x[i:(i+1)]

mat=cbind(rep(1,2),x_vec)

alpha=solve(mat)%*%u_vec

a_b[,i]=alpha

u_f=function(x){
 
N=c((x_vec[2]-x)/(x_vec[2]-x_vec[1]),(x-x_vec[1])/(x_vec[2]-x_vec[1]))

return(lam*sum(N*u_vec))
 
}
 
u_f2=function(x){

l_k=c(-1/(x_vec[2]-x_vec[1]),1/(x_vec[2]-x_vec[1]))

N=c((x_vec[2]-x)/(x_vec[2]-x_vec[1]),(x-x_vec[1])/(x_vec[2]-x_vec[1]))

return(as.numeric(l_k%*%t(t(u_vec)))*as.numeric(N%*%u_vec))

}

# monte carlo integral

Y1=0;Y2=0

n=100000

x1=sample(seq(x_vec[1],x_vec[2],0.001),n,replace=T);x2=sample(seq(x_vec[1],x_vec[2],0.001),n,replace=T)

for(j in 1:(n)){
 
Y1=Y1+u_f(x1[j]);

Y2=Y2+u_f2(x2[j]) 
 
}

# monte carlo積分値
print(Y1/n);print(Y2/n)

integral=c(integral,Y1/n+Y2/n) 
 
}


u_dif=u[2:length(u)]-u[1:(length(u)-1)]

cor(integral,u_dif)

# p.48

x=c(0,1,2);y=c(1,0,2)

pthi=sample(1:10,3)/10

mat=cbind(rep(1,3),x,y)

a=solve(mat)%*%pthi

delta=det(mat)/2

A1=(x[2]*y[3]-x[3]*y[2])/(2*delta)

A2=(x[3]*y[1]-x[1]*y[3])/(2*delta)

A3=(x[1]*y[2]-x[2]*y[1])/(2*delta)

B1=(y[2]-y[3])/(2*delta)

B2=(y[3]-y[1])/(2*delta)

B3=(y[1]-y[2])/(2*delta)

C1=(x[3]-x[2])/(2*delta)

C2=(x[1]-x[3])/(2*delta)

C3=(x[2]-x[1])/(2*delta)

# N1==1

A1+B1*x[1]+C1*y[1]

# N2==1

A2+B2*x[2]+C2*y[2]

# N3==1

A3+B3*x[3]+C3*y[3]



B=c(B1,B2,B3);C=c(C1,C2,C3)

# (3.33)

G=(t(t(B))%*%t(B)+t(t(C))%*%t(C))*delta


# 重みの最適化(pthi_e,pthi_k)でよいか？

Q_k=c(10,20)

# 境界をc(0,1)->c(2,2)

L_k=sqrt((x[3]-x[2])^2+(y[3]-y[2])^2)

H=L_k*matrix(c(2,1,1,2),ncol=2)/6

pthi_e=rep(1,3);pthi_k=rep(1,2)

cost=function(pthi1,pthi2){
 
return((as.numeric(t(c(pthi1))%*%G%*%pthi)-as.numeric(t(c(pthi2))%*%H%*%Q_k))^2) 
 
}


ite=100

r=0.001

for(l in 1:ite){
 
vec=pthi_e

for(j in 1:length(pthi_e)){
 
vec_sub=vec;vec_sub[j]=vec_sub[j]+0.01

pthi_e[j]=pthi_e[j]-r*(cost(vec_sub,pthi_k)-cost(vec,pthi_k))/0.01
 
}

vec=pthi_k

for(j in 1:length(pthi_k)){
 
vec_sub=vec;vec_sub[j]=vec_sub[j]+h

pthi_k[j]=pthi_k[j]-r*(cost(pthi_e,vec_sub)-cost(pthi_e,vec))/0.01
 
}
 
print(cost(pthi_e,pthi_k))

}