Go to Qiita Advent Calendar 2024 Top
LoginSignup
3
4

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 1 year has passed since last update.

@kozakai-ryouta(亮太 小酒井)

計算力学の基礎　丸善株式会社　八田夏夫先生

Last updated at Posted at 2021-10-03

#二重積分　p.50
#リーマン和

f=function(x,y){
return(x-y)
}

a=0
b=1
c=1
d=2

n=1000
m=1000

X=seq(a,b,1/n)
Y=seq(c,d,1/m)

value=c()
for(i in 1:length(X)){
value=c(value,f(X[i],Y))  
}

sum(value)/(n*m)



#積分領域が関数で表現される場合 p.55
#簡単な二重積分の計算例(i)

#上側
x1=function(x){
return(sqrt(1-x^2))  
}

#下側
x2=function(x){
return(0)  
}

#関数
f=function(x,y){
return(x*y)  
}

n=1000

#y軸積分領域
a=0
b=1

Y=a+(b-a)*seq(1,n,1)/n
value=0
for(i in 1:length(Y)){
B=x1(Y[i])    
A=x2(Y[i])  
if(round(B,7)!=round(A,7)){
#n;差分メッシュ
#積分の近似的計算：微分積分要論　学術図書出版　戸田暢茂先生　p.95の例4参照
seq=A+(B-A)*seq(1,n,1)/n
value2=f(seq,Y[i])
value=value+sum(value2)*(B-A)/n
}
}

#二重積分値
sum(value)*(b-a)/n




#p.65 2重積分の計算例
#変数変換を用いた例

x=function(r,sita){
return(r*cos(sita))  
}

y=function(r,sita){
return(r*sin(sita))  
}

f=function(X,Y){
return(exp(-X^2-Y^2))  
}

r=10
sita=2*pi

r_mesh=seq(0,r,0.01)
sita_mesh=seq(0,sita,1/(2*pi))
h=10^(-5)

Integral_value=0
for(i in 1:length(r_mesh)){
for(j in 1:length(sita_mesh)){  
#ヤコビアンの計算  
dxdr=(x(r_mesh[i]+h,sita_mesh[j])-x(r_mesh[i],sita_mesh[j]))/h
dxdsita=(x(r_mesh[i],sita_mesh[j]+h)-x(r_mesh[i],sita_mesh[j]))/h
dydr=(y(r_mesh[i]+h,sita_mesh[j])-y(r_mesh[i],sita_mesh[j]))/h
dydsita=(y(r_mesh[i],sita_mesh[j]+h)-y(r_mesh[i],sita_mesh[j]))/h
Jacobian=matrix(c(dxdr,dydr,dxdsita,dydsita),nrow=2)  
#積分値の合計
Integral_value=Integral_value+f(x(r_mesh[i],sita_mesh[j]),y(r_mesh[i],sita_mesh[j]))*det(Jacobian)  
}}

#数値計算による実際の積分値
sita*r*Integral_value/(length(r_mesh)*length(sita_mesh))



#p.81 3.4 曲線とその接線、法線、順法線

ax=1;bx=2
ay=2;by=3
az=2;bz=1

x=function(t){
return(ax+bx*t^2)
}

y=function(t){
return(ay+by*t)  
}

z=function(t){
return(az+bz*t^3)
}

t0=2
l=0.001

dx=(x(t0+l)-x(t0))/l
dy=(y(t0+l)-y(t0))/l
dz=(z(t0+l)-z(t0))/l

r=c(dx,dy,dz)
point=c(x(t0),y(t0),z(t0))


#Broyden Quasi-Newton Algorithm p.213

f=function(XX){
eq=(c(XX)-point)/r  
return(c(eq[1]-eq[2],eq[2]-eq[3],eq[3]-eq[1]))
}

X=rep(5,3)

ite=10^4
H=diag(f(X))
eta=10^(-3)

for(l in 1:ite){
X_pre=X  
X=X-eta*H%*%f(X)  
s=X-X_pre
y=f(X)-f(X_pre)
H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H
print(sum(abs(f(X))))
}

#接線の方程式を満たす
print((X-point)/r)

#接線を表示
xs=seq(-100,100,0.1)
ys=point[2]+(xs-point[1])*r[2]/r[1]
zs=point[3]+(xs-point[1])*r[3]/r[1]

library(rgl)
plot3d(xs,ys,zs,col=2,type="p")




#p.82 線素

ax=1;bx=2
ay=2;by=3
az=2;bz=1
l=0.001
#積分区間の決定
a=0;b=2

x=function(t){
return(ax+bx*t^2)
}

y=function(t){
return(ay+by*t)  
}

z=function(t){
return(az+bz*t^3)
}


#f(t)=sqrt(dx(t)^2+dy(t)^2+dz(t)^2)とする
f=function(t){
dx=(x(t+l)-x(t))/l
dy=(y(t+l)-y(t))/l
dz=(z(t+l)-z(t))/l
r=c(dx,dy,dz)
point=c(x(t),y(t),z(t))
return(sqrt(sum(r^2)))
}

#n;差分メッシュ
#積分の近似的計算：微分積分要論　学術図書出版　戸田暢茂先生　p.95の例4参照
n=10^5
mesh=a+(b-a)*seq(1,n,1)/n
value=c()
for(j in 1:n){
value=c(value,f(mesh[j]))
}

#線積分の結果 (p.82 (3.29))
sum(value)*(b-a)/n



#p.82 線素

#つるまき線の例

ax=1
cx=2

#積分区間の決定
a=0;b=2

#a~bまでの時刻の細分(dt)
n=10^5
dt=(b-a)/n

x=function(t){
return(ax*cos(t))
}

y=function(t){
return(ax*sin(t))  
}

z=function(t){
return(cx*t)
}


#f(t)=sqrt(dx(t)^2+dy(t)^2+dz(t)^2)とする
f=function(t){
dx=(x(t+dt)-x(t))/dt
dy=(y(t+dt)-y(t))/dt
dz=(z(t+dt)-z(t))/dt
r=c(dx,dy,dz)
point=c(x(t),y(t),z(t))
return(sqrt(sum(r^2)))
}

#n;差分メッシュ
#積分の近似的計算：微分積分要論　学術図書出版　戸田暢茂先生　p.95の例4参照
mesh=a+(b-a)*seq(1,n,1)/n
value=c()
for(j in 1:n){
value=c(value,f(mesh[j]))
}

#線積分の結果 (p.82 (3.29))
sum(value)*(b-a)/n

#積分の値
sqrt(ax^2+cx^2)*(b-a)
 
t=2
dx=(x(t+dt)-x(t))/dt
dy=(y(t+dt)-y(t))/dt
dz=(z(t+dt)-z(t))/dt
r=c(dx,dy,dz)
tt=c(dx,dy,dz)/sqrt(sum(r^2))

#ttと同じになるはず
c(-ax*sin(t),ax*cos(t),cx)/sqrt(sum(r^2))

d2x=(x(t+dt)-2*x(t)+x(t-dt))/(dt^2)
d2y=(y(t+dt)-2*y(t)+y(t-dt))/(dt^2)
d2z=(z(t+dt)-2*z(t)+z(t-dt))/(dt^2)
tdash=c(d2x,d2y,d2z)/sqrt(sum(r^2))

#tdashと同じになるはず
c(-ax*cos(t),-ax*sin(t),0)/sqrt(sum(r^2))

#長さを1にする
n=tdash/sqrt(sum(tdash^2))

mat=as.matrix(rbind(r,n))
bb=c(det(mat[,2:3]),det(mat[,c(3,1)]),det(mat[,1:2]))/sqrt(sum(r^2))

#bbと同じになるはず
c(cx*sin(t),-cx*cos(t),ax)/sqrt(sum(r^2))



#3.5 曲面と曲面積 p.86

aa=1

x=function(u,v){
return(aa*sin(u)*cos(v))  
}

y=function(u,v){
return(aa*sin(u)*sin(v))  
}

z=function(u,v){
return(aa*cos(u))  
}

#上側
x1=function(s){
return(pi)  
}

#下側
x2=function(s){
return(0)  
}

#関数
h=0.001
f=function(u,v){
r1=c((x(u+h,v)-x(u,v))/h,(y(u+h,v)-y(u,v))/h,(z(u+h,v)-z(u,v))/h)
r2=c((x(u,v+h)-x(u,v))/h,(y(u,v+h)-y(u,v))/h,(z(u,v+h)-z(u,v))/h)
e=sum(r1*r1)
g=sum(r2*r2)
f=sum(r1*r2)
return(sqrt(e*g-f^2))  
}

n=1000

#y軸積分領域
a=0
b=2*pi

Y=a+(b-a)*seq(1,n,1)/n
value=c()
for(i in 1:length(Y)){
B=x1(Y[i])    
A=x2(Y[i])  
if(round(B,7)!=round(A,7)){
#n;差分メッシュ
#積分の近似的計算：微分積分要論　学術図書出版　戸田暢茂先生　p.95の例4参照
seq=A+(B-A)*seq(1,n,1)/n
value2=c()
for(j in 1:length(seq)){
value2=c(value2,f(seq[j],Y[i]))}
value=c(value,sum(value2)*(B-A)/n)
}
}

#二重積分値
sum(value)*(b-a)/n



#オイラー法　p.118

#空気中の水滴の落下運動(抗力を考慮した場合)

#落下速度(m/s)
v=0.001
#分割時間(s)
tau=0.01
#水滴半径の増加率(m/s)
a=10^(-7)
#水滴半径(μ/m)
r=3
#水の密度(kg/m3)
pL=1000
#空気密度(kg/m3)
pa=1.2
#重力加速度(m/s2)
g=9.8
#空気の動粘性係数(m2/s)
nyu=1.5*10^(-5)
#レイノルズ数
Re=2*r*v/nyu
#水滴を球体とした場合の抗力係数
C=24*(1+0.15*Re^(0.687))/Re

ite=20000
v_vec=c()
r_vec=c()

for(i in 1:ite){
v=v+(g/a-(3/r)*(v+C*pa*v^2/(8*a*pL)))*a*tau  
r=r+a*tau  
v_vec=c(v_vec,v)
r_vec=c(r_vec,r)
}

#水滴の半径(r)の上昇により、落下速度が減少、一定に落ち着く(重力と抵抗力がつり合って、水滴の終端速度に落ち着く)
plot(r_vec,v_vec)



#p.127 オイラー法

#空気中の水滴の落下運動で、x軸方向から一定の風速で流れている場合(抗力を考慮した場合)

#落下速度(m/s)
v=0.001
#分割時間(s)
tau=0.01
#水滴半径の増加率(m/s)
a=10^(-7)
#水滴半径(μ/m)
r=3
#水の密度(kg/m3)
pL=1000
#空気密度(kg/m3)
pa=1.2
#重力加速度(m/s2)
g=9.8
#空気の動粘性係数(m2/s)
nyu=1.5*10^(-5)

#水滴を球体とした場合の抗力係数
C=24*(1+0.15*Re^(0.687))/Re
#流体の流れている速度(U)
U=2
#x軸初期位置
x=0
#y軸初期位置
y=0

dx=0;dy=0

#反復回数
ite=10000

x_vec=c()
y_vec=c()
r_vec=c()
v_vec=c()

for(l in 1:ite){
  
#合相対速度
v=sqrt((U-dx)^2+dy^2)

#レイノルズ数
Re=2*r*v/nyu

#変化量
dx=dx-3*a*dx/r+(3*C*pa*sqrt((U-dx)^2+dy^2)*(U-dx))/(8*r*pL)*tau  
dy=dy-3*a*dy/r+g-(3*C*pa*sqrt((U-dx)^2+dy^2)*(dy))/(8*r*pL)*tau

#座標、水滴半径の更新
x=x+dx*tau
y=y+dy*tau
r=r+a*tau

#データの保存
x_vec=c(x_vec,x)
y_vec=c(y_vec,y)
r_vec=c(r_vec,r)
v_vec=c(v_vec,v)

}

#水滴の半径(r)の上昇により、落下速度が減少、一定に落ち着く(重力と抵抗力がつり合って、水滴の終端速度に落ち着く)
plot(r_vec,v_vec)



#p.185 曲面の曲率の極小値、極大値及び両者の和

z=function(x,y){
return(3*x^2+y^2+x*y)  
}


X=20;Y=40
l=0.001

p=(z(X+l,Y)-z(X,Y))/l
q=(z(X,Y+l)-z(X,Y))/l
r=(z(X+2*l,Y)-2*z(X+l,Y)+z(X,Y))/(l^2)
h=(z(X+l,Y+l)-z(X+l,Y)-z(X,Y+l)+z(X,Y))/(l^2)
t=(z(X,Y+2*l)-2*z(X,Y+l)+z(X,Y))/(l^2)
rx=c(1,0,p)
ry=c(0,1,q)

a=-p
b=-q
c=1

e=1+p^2
f=p*q
g=1+q^2

#外積の大きさ
e*g-f^2

A=a/sqrt(e*g-f^2)
B=b/sqrt(e*g-f^2)
C=c/sqrt(e*g-f^2)
L=r/sqrt(e*g-f^2)
M=h/sqrt(e*g-f^2)
N=t/sqrt(e*g-f^2)

#曲率半径の極値をとる角度
sita1=atan(2*h/(r-t))
sita2=sita1+pi/2

#極大、極小値の曲率半径
R1=1/((r+t)/2+(r-t)*cos(2*sita1)/2+h*sin(2*sita1))
R2=1/((r+t)/2+(r-t)*cos(2*sita2)/2+h*sin(2*sita2))

#オイラーの公式
Euler=function(sita){
return(cos(sita)^2/R1+sin(sita)^2/R2)  
}

Euler(0)　#確認1
Euler(pi/2)　#確認2

#直交する2つの法平面によって切断したとき、2つの切断面に現れる曲率の曲率の和は不変
s=1.5　#任意のsについて成立
Euler(s)+Euler(s+pi/2)==1/R1+1/R2

#主方向を解の公式で計算　p.191
aa=f*N-g*M
bb=e*N-g*L
cc=e*M-f*L
U1=(-bb+sqrt(bb^2-4*aa*cc))/(2*aa)
U2=(-bb-sqrt(bb^2-4*aa*cc))/(2*aa)

#解になっているかどうか確認
aa*U1^2+bb*U1+cc
aa*U2^2+bb*U2+cc




#p.192

#f(x,y,z)=uについて

#直交しない場合

z=function(x,y){
return(x^2+y^2-3*x*y)  
}

u=function(x,y,l){
return(x^2+y^2+l^2)
}

v=function(x,y,l){
return(x^2+y^2-l^3)  
}

w=function(x,y,l){
return(x^2+y^2+l^3)  
}


h=0.001
X=1;Y=5

dux=(u(X+h,Y,z(X,Y))-u(X,Y,z(X,Y)))/h
duy=(u(X,Y+h,z(X,Y))-u(X,Y,z(X,Y)))/h
duz=(u(X,Y,z(X,Y)+h)-u(X,Y,z(X,Y)))/h

dvx=(v(X+h,Y,z(X,Y))-v(X,Y,z(X,Y)))/h
dvy=(v(X,Y+h,z(X,Y))-v(X,Y,z(X,Y)))/h
dvz=(v(X,Y,z(X,Y)+h)-v(X,Y,z(X,Y)))/h

dwx=(w(X+h,Y,z(X,Y))-w(X,Y,z(X,Y)))/h
dwy=(w(X,Y+h,z(X,Y))-w(X,Y,z(X,Y)))/h
dwz=(w(X,Y,z(X,Y)+h)-w(X,Y,z(X,Y)))/h


h1=sqrt(dux^2+duy^2+duz^2)
h2=sqrt(dvx^2+dvy^2+dvz^2)
h3=sqrt(dwx^2+dwy^2+dwz^2)

lam=c(dux,dvx,dwx)/c(h1,h2,h3)
mu=c(duy,dvy,dwy)/c(h1,h2,h3)
nyu=c(duz,dvz,dwz)/c(h1,h2,h3)

#直交条件
lam[1]*lam[2]+mu[1]*mu[2]+nyu[1]*nyu[2]
lam[2]*lam[3]+mu[2]*mu[3]+nyu[2]*nyu[3]
lam[1]*lam[3]+mu[1]*mu[3]+nyu[1]*nyu[3]




#p.192

#f(x,y,z)=uについて

#直交する場合

z=function(x,y){
return(x^2+y^2-3*x*y)  
}

u=function(x,y,l){
return(l^2)
}

v=function(x,y,l){
return(x^2)  
}

w=function(x,y,l){
return(y^2)  
}


h=0.001
X=1;Y=5

dux=(u(X+h,Y,z(X,Y))-u(X,Y,z(X,Y)))/h
duy=(u(X,Y+h,z(X,Y))-u(X,Y,z(X,Y)))/h
duz=(u(X,Y,z(X,Y)+h)-u(X,Y,z(X,Y)))/h

dvx=(v(X+h,Y,z(X,Y))-v(X,Y,z(X,Y)))/h
dvy=(v(X,Y+h,z(X,Y))-v(X,Y,z(X,Y)))/h
dvz=(v(X,Y,z(X,Y)+h)-v(X,Y,z(X,Y)))/h

dwx=(w(X+h,Y,z(X,Y))-w(X,Y,z(X,Y)))/h
dwy=(w(X,Y+h,z(X,Y))-w(X,Y,z(X,Y)))/h
dwz=(w(X,Y,z(X,Y)+h)-w(X,Y,z(X,Y)))/h


h1=sqrt(dux^2+duy^2+duz^2)
h2=sqrt(dvx^2+dvy^2+dvz^2)
h3=sqrt(dwx^2+dwy^2+dwz^2)

lam=c(dux,dvx,dwx)/c(h1,h2,h3)
mu=c(duy,dvy,dwy)/c(h1,h2,h3)
nyu=c(duz,dvz,dwz)/c(h1,h2,h3)

#直交条件
lam[1]*lam[2]+mu[1]*mu[2]+nyu[1]*nyu[2]
lam[2]*lam[3]+mu[2]*mu[3]+nyu[2]*nyu[3]
lam[1]*lam[3]+mu[1]*mu[3]+nyu[1]*nyu[3]




#極座標での計算 p.199

x=function(r,s,l){
return(r*sin(s)*cos(l)) 
}

y=function(r,s,l){
return(r*sin(s)*sin(l))  
}

z=function(r,s,l){
return(r*cos(s))  
}

R=3
sita=pi/4
pthi=pi/3
h=0.0001

dxr=(x(R+h,sita,pthi)-x(R,sita,pthi))/h
dxs=(x(R,sita+h,pthi)-x(R,sita,pthi))/h
dxl=(x(R,sita,pthi+h)-x(R,sita,pthi))/h

dyr=(y(R+h,sita,pthi)-y(R,sita,pthi))/h
dys=(y(R,sita+h,pthi)-y(R,sita,pthi))/h
dyl=(y(R,sita,pthi+h)-y(R,sita,pthi))/h

dzr=(z(R+h,sita,pthi)-z(R,sita,pthi))/h
dzs=(z(R,sita+h,pthi)-z(R,sita,pthi))/h
dzl=(z(R,sita,pthi+h)-z(R,sita,pthi))/h

g1=sqrt(dxr^2+dyr^2+dzr^2)
g2=sqrt(dxs^2+dys^2+dzs^2)
g3=sqrt(dxl^2+dyl^2+dzl^2)

1             #g1
R             #g2
R*sin(sita)   #g3



#p.201
#直交座標系から平面極座標系へ

#6.1 直交曲線座標系における基本単位ベクトルの考え方
r=3
sita=2*pi/3
h=0.0001

guzai1=r
guzai2=sita

x1=function(R,SITA){
return(R*cos(SITA))
}

x2=function(R,SITA){
return(R*sin(SITA))  
}

dx1_dguzai1=(x1(guzai1+h,guzai2)-x1(guzai1,guzai2))/h
dx2_dguzai1=(x2(guzai1+h,guzai2)-x2(guzai1,guzai2))/h
dx1_dguzai2=(x1(guzai1,guzai2+h)-x1(guzai1,guzai2))/h
dx2_dguzai2=(x2(guzai1,guzai2+h)-x2(guzai1,guzai2))/h

g1=sqrt(dx1_dguzai1^2+dx2_dguzai1^2)
g2=sqrt(dx1_dguzai2^2+dx2_dguzai2^2)

#直交座標系での単位ベクトル
e1hat=c(1,0)
e2hat=c(0,1)

#平面極座標系での単位ベクトル
e1=(dx1_dguzai1*e1hat+dx2_dguzai1*e2hat)/g1
e2=(dx1_dguzai2*e1hat+dx2_dguzai2*e2hat)/g2




#p.205
#6.2 直交曲線座標系における勾配、発散、回転などの一般的表示

#6.3 円柱座標

x1=function(X,R,SITA){
return(X)  
}

x2=function(X,R,SITA){
return(R*cos(SITA))  
}

x3=function(X,R,SITA){
return(R*sin(SITA))  
}

h=0.001
x=2
r=3
sita=2*pi/3

e_g=function(X,R,SITA){

guzai1=X
guzai2=R
guzai3=SITA

dx1_dguzai1=(x1(guzai1+h,guzai2,guzai3)-x1(guzai1,guzai2,guzai3))/h
dx1_dguzai2=(x1(guzai1,guzai2+h,guzai3)-x1(guzai1,guzai2,guzai3))/h
dx1_dguzai3=(x1(guzai1,guzai2,guzai3+h)-x1(guzai1,guzai2,guzai3))/h

dx2_dguzai1=(x2(guzai1+h,guzai2,guzai3)-x2(guzai1,guzai2,guzai3))/h
dx2_dguzai2=(x2(guzai1,guzai2+h,guzai3)-x2(guzai1,guzai2,guzai3))/h
dx2_dguzai3=(x2(guzai1,guzai2,guzai3+h)-x2(guzai1,guzai2,guzai3))/h

dx3_dguzai1=(x3(guzai1+h,guzai2,guzai3)-x3(guzai1,guzai2,guzai3))/h
dx3_dguzai2=(x3(guzai1,guzai2+h,guzai3)-x3(guzai1,guzai2,guzai3))/h
dx3_dguzai3=(x3(guzai1,guzai2,guzai3+h)-x3(guzai1,guzai2,guzai3))/h

g1=sqrt(dx1_dguzai1^2+dx2_dguzai1^2+dx3_dguzai1^2)
g2=sqrt(dx1_dguzai2^2+dx2_dguzai2^2+dx3_dguzai2^2)
g3=sqrt(dx1_dguzai3^2+dx2_dguzai3^2+dx3_dguzai3^2)

#直交座標系での単位ベクトル
e1hat=c(1,0,0)
e2hat=c(0,1,0)
e3hat=c(0,0,1)

#直交曲線座標系での単位ベクトル
e1=(dx1_dguzai1*e1hat+dx2_dguzai1*e2hat+dx3_dguzai1*e3hat)/g1
e2=(dx1_dguzai2*e1hat+dx2_dguzai2*e2hat+dx3_dguzai2*e3hat)/g2
e3=(dx1_dguzai3*e1hat+dx2_dguzai3*e2hat+dx3_dguzai3*e3hat)/g3

return(cbind(e1,e2,e3,c(g1,g2,g3)))
}

A=function(X,R,SITA){
return(c(X^2,R,sin(SITA)))
}

guzai1=x
guzai2=r
guzai3=sita

Avec=A(guzai1,guzai2,guzai3)
dA_dguzai1=(A(guzai1+h,guzai2,guzai3)-A(guzai1,guzai2,guzai3))/h
dA_dguzai2=(A(guzai1,guzai2+h,guzai3)-A(guzai1,guzai2,guzai3))/h
dA_dguzai3=(A(guzai1,guzai2,guzai3+h)-A(guzai1,guzai2,guzai3))/h

e_g_mat=e_g(guzai1,guzai2,guzai3)
e1=e_g_mat[,1]
e2=e_g_mat[,2]
e3=e_g_mat[,3]
g=e_g_mat[,4]
g1=g[1]
g2=g[2]
g3=g[3]

e_g_dguzai1=e_g(guzai1+h,guzai2,guzai3)
e_g_dguzai2=e_g(guzai1,guzai2+h,guzai3)
e_g_dguzai3=e_g(guzai1,guzai2,guzai3+h)

g_dguzai1=e_g_dguzai1[,4]
g_dguzai2=e_g_dguzai2[,4]
g_dguzai3=e_g_dguzai3[,4]


#勾配
grad_A=e1*dA_dguzai1/g1+e2*dA_dguzai2/g2+e3*dA_dguzai3/g3

#回転
Wx1=((g_dguzai2[3]*A(guzai1,guzai2+h,guzai3)[3]-g3*Avec[3])/h-(g_dguzai3[2]*A(guzai1,guzai2,guzai3+h)[2]-g2*Avec[2])/h)/(g2*g3)

Wx2=((g_dguzai3[1]*A(guzai1,guzai2,guzai3+h)[1]-g1*Avec[1])/h-(g_dguzai1[3]*A(guzai1+h,guzai2,guzai3)[3]-g3*Avec[3])/h)/(g1*g3)

Wx3=((g_dguzai1[2]*A(guzai1+h,guzai2,guzai3)[2]-g2*Avec[2])/h-(g_dguzai2[1]*A(guzai1,guzai2+h,guzai3)[1]-g1*Avec[1])/h)/(g2*g1)

rot_A=e1*Wx1+e2*Wx2+e3*Wx3



#p.213
#6.2 直交曲線座標系における勾配、発散、回転などの一般的表示

#6.3.3 球面座標

x1=function(R,SITA,X){
return(R*cos(SITA))  
}

x2=function(R,SITA,X){
return(R*sin(SITA)*cos(X))  
}

x3=function(R,SITA,X){
return(R*sin(X)*sin(SITA))  
}

h=0.00001
x=pi/3
r=3
sita=2*pi/3

e_g=function(R,SITA,X){

guzai1=R
guzai2=SITA
guzai3=X

dx1_dguzai1=(x1(guzai1+h,guzai2,guzai3)-x1(guzai1,guzai2,guzai3))/h
dx1_dguzai2=(x1(guzai1,guzai2+h,guzai3)-x1(guzai1,guzai2,guzai3))/h
dx1_dguzai3=(x1(guzai1,guzai2,guzai3+h)-x1(guzai1,guzai2,guzai3))/h

dx2_dguzai1=(x2(guzai1+h,guzai2,guzai3)-x2(guzai1,guzai2,guzai3))/h
dx2_dguzai2=(x2(guzai1,guzai2+h,guzai3)-x2(guzai1,guzai2,guzai3))/h
dx2_dguzai3=(x2(guzai1,guzai2,guzai3+h)-x2(guzai1,guzai2,guzai3))/h

dx3_dguzai1=(x3(guzai1+h,guzai2,guzai3)-x3(guzai1,guzai2,guzai3))/h
dx3_dguzai2=(x3(guzai1,guzai2+h,guzai3)-x3(guzai1,guzai2,guzai3))/h
dx3_dguzai3=(x3(guzai1,guzai2,guzai3+h)-x3(guzai1,guzai2,guzai3))/h

g1=sqrt(dx1_dguzai1^2+dx2_dguzai1^2+dx3_dguzai1^2)
g2=sqrt(dx1_dguzai2^2+dx2_dguzai2^2+dx3_dguzai2^2)
g3=sqrt(dx1_dguzai3^2+dx2_dguzai3^2+dx3_dguzai3^2)

#直交座標系での単位ベクトル
e1hat=c(1,0,0)
e2hat=c(0,1,0)
e3hat=c(0,0,1)

#直交曲線座標系での単位ベクトル
e1=(dx1_dguzai1*e1hat+dx2_dguzai1*e2hat+dx3_dguzai1*e3hat)/g1
e2=(dx1_dguzai2*e1hat+dx2_dguzai2*e2hat+dx3_dguzai2*e3hat)/g2
e3=(dx1_dguzai3*e1hat+dx2_dguzai3*e2hat+dx3_dguzai3*e3hat)/g3

return(cbind(e1,e2,e3,c(g1,g2,g3)))
}

A=function(R,SITA,X){
return(c(X^2,R,sin(SITA)))
}

guzai1=r
guzai2=sita
guzai3=x

Avec=A(guzai1,guzai2,guzai3)
dA_dguzai1=(A(guzai1+h,guzai2,guzai3)-A(guzai1,guzai2,guzai3))/h
dA_dguzai2=(A(guzai1,guzai2+h,guzai3)-A(guzai1,guzai2,guzai3))/h
dA_dguzai3=(A(guzai1,guzai2,guzai3+h)-A(guzai1,guzai2,guzai3))/h

e_g_mat=e_g(guzai1,guzai2,guzai3)
e1=e_g_mat[,1]
e2=e_g_mat[,2]
e3=e_g_mat[,3]
g=e_g_mat[,4]
g1=g[1]
g2=g[2]
g3=g[3]

e_g_dguzai1=e_g(guzai1+h,guzai2,guzai3)
e_g_dguzai2=e_g(guzai1,guzai2+h,guzai3)
e_g_dguzai3=e_g(guzai1,guzai2,guzai3+h)

g_dguzai1=e_g_dguzai1[,4]
g_dguzai2=e_g_dguzai2[,4]
g_dguzai3=e_g_dguzai3[,4]


#勾配
grad_A=e1*dA_dguzai1/g1+e2*dA_dguzai2/g2+e3*dA_dguzai3/g3

#回転
Wx1=((g_dguzai2[3]*A(guzai1,guzai2+h,guzai3)[3]-g3*Avec[3])/h-(g_dguzai3[2]*A(guzai1,guzai2,guzai3+h)[2]-g2*Avec[2])/h)/(g2*g3)

Wx2=((g_dguzai3[1]*A(guzai1,guzai2,guzai3+h)[1]-g1*Avec[1])/h-(g_dguzai1[3]*A(guzai1+h,guzai2,guzai3)[3]-g3*Avec[3])/h)/(g1*g3)

Wx3=((g_dguzai1[2]*A(guzai1+h,guzai2,guzai3)[2]-g2*Avec[2])/h-(g_dguzai2[1]*A(guzai1,guzai2+h,guzai3)[1]-g1*Avec[1])/h)/(g2*g1)

rot_A=e1*Wx1+e2*Wx2+e3*Wx3



#p.214 ひずみとひずみ速度

u=function(x,i){
if(i==1){  
return(x[1]+x[2]^2+x[3])  
}else{
if(i==2){
return(x[1]+x[2]+x[3])  
}else{
return(x[1]^2+x[2]^2+x[3])  
}}}



ep=array(0,dim=c(3,3))
abc=ep

h=0.001

X=c(4,3,6)

for(k in 1:ncol(ep)){
for(j in 1:nrow(ep)){
if(k==j){ 
vector=X;vector_h=vector
vector_h[j]=vector_h[j]+h
ep[j,k]=(u(vector_h,k)-u(vector,k))/h  
}else{
vector1=X;vector_h1=vector1
vector_h1[j]=vector_h1[j]+h
vector2=X;vector_h2=vector2
vector_h2[k]=vector_h2[k]+h
ep[j,k]=((u(vector_h1,k)-u(vector1,k))/h+(u(vector_h2,j)-u(vector2,j))/h)/2
abc[j,k]=((u(vector_h1,k)-u(vector1,k))/h-(u(vector_h2,j)-u(vector2,j))/h)/2
}}}

dX=c(0.01,0.2,0.07)

alpha=abc[3,1]
beta=abc[1,2]
gamma=abc[2,3]

du=ep%*%X

du[1]=du[1]+(alpha*dX[3]-beta*dX[2])
du[2]=du[2]+(beta*dX[1]-gamma*dX[3])
du[3]=du[3]+(gamma*dX[2]-alpha*dX[1])
3
4
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
3
4

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?