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

```r

#二重積分　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)

```

```r

#積分領域が関数で表現される場合 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


```

```r

#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))

```


```r

#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")


```

```r

#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

```

```r

#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))

```

```r

#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

```

```r

#オイラー法　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)

```


```r

#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)

```



```r

#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


```

```r

#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]


```


```r

#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]


```


```r

#極座標での計算 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

```

```r

#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


```

```r

#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

```

```r

#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

```


```r

#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])


```