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

#二重積分　p.50#リーマン和f=function(x,y){return(x-y)}a=0b=1c=1d=2n=1000m=1000X=seq(a,b,1/n)Y=seq(c,d,1/m)value=c()for(iin1: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=0b=1Y=a+(b-a)*seq(1,n,1)/nvalue=0for(iin1: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)/nvalue2=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=10sita=2*pir_mesh=seq(0,r,0.01)sita_mesh=seq(0,sita,1/(2*pi))h=10^(-5)Integral_value=0for(iin1:length(r_mesh)){for(jin1:length(sita_mesh)){#ヤコビアンの計算 dxdr=(x(r_mesh[i]+h,sita_mesh[j])-x(r_mesh[i],sita_mesh[j]))/hdxdsita=(x(r_mesh[i],sita_mesh[j]+h)-x(r_mesh[i],sita_mesh[j]))/hdydr=(y(r_mesh[i]+h,sita_mesh[j])-y(r_mesh[i],sita_mesh[j]))/hdydsita=(y(r_mesh[i],sita_mesh[j]+h)-y(r_mesh[i],sita_mesh[j]))/hJacobian=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=2ay=2;by=3az=2;bz=1x=function(t){return(ax+bx*t^2)}y=function(t){return(ay+by*t)}z=function(t){return(az+bz*t^3)}t0=2l=0.001dx=(x(t0+l)-x(t0))/ldy=(y(t0+l)-y(t0))/ldz=(z(t0+l)-z(t0))/lr=c(dx,dy,dz)point=c(x(t0),y(t0),z(t0))#Broyden Quasi-Newton Algorithm p.213f=function(XX){eq=(c(XX)-point)/rreturn(c(eq[1]-eq[2],eq[2]-eq[3],eq[3]-eq[1]))}X=rep(5,3)ite=10^4H=diag(f(X))eta=10^(-3)for(lin1:ite){X_pre=XX=X-eta*H%*%f(X)s=X-X_prey=f(X)-f(X_pre)H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%Hprint(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=2ay=2;by=3az=2;bz=1l=0.001#積分区間の決定a=0;b=2x=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))/ldy=(y(t+l)-y(t))/ldz=(z(t+l)-z(t))/lr=c(dx,dy,dz)point=c(x(t),y(t),z(t))return(sqrt(sum(r^2)))}#n;差分メッシュ#積分の近似的計算：微分積分要論　学術図書出版　戸田暢茂先生　p.95の例4参照n=10^5mesh=a+(b-a)*seq(1,n,1)/nvalue=c()for(jin1:n){value=c(value,f(mesh[j]))}#線積分の結果 (p.82 (3.29))sum(value)*(b-a)/n

#p.82 線素#つるまき線の例ax=1cx=2#積分区間の決定a=0;b=2#a~bまでの時刻の細分(dt)n=10^5dt=(b-a)/nx=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))/dtdy=(y(t+dt)-y(t))/dtdz=(z(t+dt)-z(t))/dtr=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)/nvalue=c()for(jin1:n){value=c(value,f(mesh[j]))}#線積分の結果 (p.82 (3.29))sum(value)*(b-a)/n#積分の値sqrt(ax^2+cx^2)*(b-a)t=2dx=(x(t+dt)-x(t))/dtdy=(y(t+dt)-y(t))/dtdz=(z(t+dt)-z(t))/dtr=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.86aa=1x=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.001f=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=0b=2*piY=a+(b-a)*seq(1,n,1)/nvalue=c()for(iin1: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)/nvalue2=c()for(jin1: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))/Reite=20000v_vec=c()r_vec=c()for(iin1:ite){v=v+(g/a-(3/r)*(v+C*pa*v^2/(8*a*pL)))*a*taur=r+a*tauv_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=0dx=0;dy=0#反復回数ite=10000x_vec=c()y_vec=c()r_vec=c()v_vec=c()for(lin1: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)*taudy=dy-3*a*dy/r+g-(3*C*pa*sqrt((U-dx)^2+dy^2)*(dy))/(8*r*pL)*tau#座標、水滴半径の更新x=x+dx*tauy=y+dy*taur=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=40l=0.001p=(z(X+l,Y)-z(X,Y))/lq=(z(X,Y+l)-z(X,Y))/lr=(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=-pb=-qc=1e=1+p^2f=p*qg=1+q^2#外積の大きさe*g-f^2A=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)#確認1Euler(pi/2)#確認2#直交する2つの法平面によって切断したとき、2つの切断面に現れる曲率の曲率の和は不変s=1.5#任意のsについて成立Euler(s)+Euler(s+pi/2)==1/R1+1/R2#主方向を解の公式で計算　p.191aa=f*N-g*Mbb=e*N-g*Lcc=e*M-f*LU1=(-bb+sqrt(bb^2-4*aa*cc))/(2*aa)U2=(-bb-sqrt(bb^2-4*aa*cc))/(2*aa)#解になっているかどうか確認aa*U1^2+bb*U1+ccaa*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.001X=1;Y=5dux=(u(X+h,Y,z(X,Y))-u(X,Y,z(X,Y)))/hduy=(u(X,Y+h,z(X,Y))-u(X,Y,z(X,Y)))/hduz=(u(X,Y,z(X,Y)+h)-u(X,Y,z(X,Y)))/hdvx=(v(X+h,Y,z(X,Y))-v(X,Y,z(X,Y)))/hdvy=(v(X,Y+h,z(X,Y))-v(X,Y,z(X,Y)))/hdvz=(v(X,Y,z(X,Y)+h)-v(X,Y,z(X,Y)))/hdwx=(w(X+h,Y,z(X,Y))-w(X,Y,z(X,Y)))/hdwy=(w(X,Y+h,z(X,Y))-w(X,Y,z(X,Y)))/hdwz=(w(X,Y,z(X,Y)+h)-w(X,Y,z(X,Y)))/hh1=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.001X=1;Y=5dux=(u(X+h,Y,z(X,Y))-u(X,Y,z(X,Y)))/hduy=(u(X,Y+h,z(X,Y))-u(X,Y,z(X,Y)))/hduz=(u(X,Y,z(X,Y)+h)-u(X,Y,z(X,Y)))/hdvx=(v(X+h,Y,z(X,Y))-v(X,Y,z(X,Y)))/hdvy=(v(X,Y+h,z(X,Y))-v(X,Y,z(X,Y)))/hdvz=(v(X,Y,z(X,Y)+h)-v(X,Y,z(X,Y)))/hdwx=(w(X+h,Y,z(X,Y))-w(X,Y,z(X,Y)))/hdwy=(w(X,Y+h,z(X,Y))-w(X,Y,z(X,Y)))/hdwz=(w(X,Y,z(X,Y)+h)-w(X,Y,z(X,Y)))/hh1=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.199x=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=3sita=pi/4pthi=pi/3h=0.0001dxr=(x(R+h,sita,pthi)-x(R,sita,pthi))/hdxs=(x(R,sita+h,pthi)-x(R,sita,pthi))/hdxl=(x(R,sita,pthi+h)-x(R,sita,pthi))/hdyr=(y(R+h,sita,pthi)-y(R,sita,pthi))/hdys=(y(R,sita+h,pthi)-y(R,sita,pthi))/hdyl=(y(R,sita,pthi+h)-y(R,sita,pthi))/hdzr=(z(R+h,sita,pthi)-z(R,sita,pthi))/hdzs=(z(R,sita+h,pthi)-z(R,sita,pthi))/hdzl=(z(R,sita,pthi+h)-z(R,sita,pthi))/hg1=sqrt(dxr^2+dyr^2+dzr^2)g2=sqrt(dxs^2+dys^2+dzs^2)g3=sqrt(dxl^2+dyl^2+dzl^2)1#g1R#g2R*sin(sita)#g3

#p.201#直交座標系から平面極座標系へ#6.1 直交曲線座標系における基本単位ベクトルの考え方r=3sita=2*pi/3h=0.0001guzai1=rguzai2=sitax1=function(R,SITA){return(R*cos(SITA))}x2=function(R,SITA){return(R*sin(SITA))}dx1_dguzai1=(x1(guzai1+h,guzai2)-x1(guzai1,guzai2))/hdx2_dguzai1=(x2(guzai1+h,guzai2)-x2(guzai1,guzai2))/hdx1_dguzai2=(x1(guzai1,guzai2+h)-x1(guzai1,guzai2))/hdx2_dguzai2=(x2(guzai1,guzai2+h)-x2(guzai1,guzai2))/hg1=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)/g1e2=(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.001x=2r=3sita=2*pi/3e_g=function(X,R,SITA){guzai1=Xguzai2=Rguzai3=SITAdx1_dguzai1=(x1(guzai1+h,guzai2,guzai3)-x1(guzai1,guzai2,guzai3))/hdx1_dguzai2=(x1(guzai1,guzai2+h,guzai3)-x1(guzai1,guzai2,guzai3))/hdx1_dguzai3=(x1(guzai1,guzai2,guzai3+h)-x1(guzai1,guzai2,guzai3))/hdx2_dguzai1=(x2(guzai1+h,guzai2,guzai3)-x2(guzai1,guzai2,guzai3))/hdx2_dguzai2=(x2(guzai1,guzai2+h,guzai3)-x2(guzai1,guzai2,guzai3))/hdx2_dguzai3=(x2(guzai1,guzai2,guzai3+h)-x2(guzai1,guzai2,guzai3))/hdx3_dguzai1=(x3(guzai1+h,guzai2,guzai3)-x3(guzai1,guzai2,guzai3))/hdx3_dguzai2=(x3(guzai1,guzai2+h,guzai3)-x3(guzai1,guzai2,guzai3))/hdx3_dguzai3=(x3(guzai1,guzai2,guzai3+h)-x3(guzai1,guzai2,guzai3))/hg1=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)/g1e2=(dx1_dguzai2*e1hat+dx2_dguzai2*e2hat+dx3_dguzai2*e3hat)/g2e3=(dx1_dguzai3*e1hat+dx2_dguzai3*e2hat+dx3_dguzai3*e3hat)/g3return(cbind(e1,e2,e3,c(g1,g2,g3)))}A=function(X,R,SITA){return(c(X^2,R,sin(SITA)))}guzai1=xguzai2=rguzai3=sitaAvec=A(guzai1,guzai2,guzai3)dA_dguzai1=(A(guzai1+h,guzai2,guzai3)-A(guzai1,guzai2,guzai3))/hdA_dguzai2=(A(guzai1,guzai2+h,guzai3)-A(guzai1,guzai2,guzai3))/hdA_dguzai3=(A(guzai1,guzai2,guzai3+h)-A(guzai1,guzai2,guzai3))/he_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.00001x=pi/3r=3sita=2*pi/3e_g=function(R,SITA,X){guzai1=Rguzai2=SITAguzai3=Xdx1_dguzai1=(x1(guzai1+h,guzai2,guzai3)-x1(guzai1,guzai2,guzai3))/hdx1_dguzai2=(x1(guzai1,guzai2+h,guzai3)-x1(guzai1,guzai2,guzai3))/hdx1_dguzai3=(x1(guzai1,guzai2,guzai3+h)-x1(guzai1,guzai2,guzai3))/hdx2_dguzai1=(x2(guzai1+h,guzai2,guzai3)-x2(guzai1,guzai2,guzai3))/hdx2_dguzai2=(x2(guzai1,guzai2+h,guzai3)-x2(guzai1,guzai2,guzai3))/hdx2_dguzai3=(x2(guzai1,guzai2,guzai3+h)-x2(guzai1,guzai2,guzai3))/hdx3_dguzai1=(x3(guzai1+h,guzai2,guzai3)-x3(guzai1,guzai2,guzai3))/hdx3_dguzai2=(x3(guzai1,guzai2+h,guzai3)-x3(guzai1,guzai2,guzai3))/hdx3_dguzai3=(x3(guzai1,guzai2,guzai3+h)-x3(guzai1,guzai2,guzai3))/hg1=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)/g1e2=(dx1_dguzai2*e1hat+dx2_dguzai2*e2hat+dx3_dguzai2*e3hat)/g2e3=(dx1_dguzai3*e1hat+dx2_dguzai3*e2hat+dx3_dguzai3*e3hat)/g3return(cbind(e1,e2,e3,c(g1,g2,g3)))}A=function(R,SITA,X){return(c(X^2,R,sin(SITA)))}guzai1=rguzai2=sitaguzai3=xAvec=A(guzai1,guzai2,guzai3)dA_dguzai1=(A(guzai1+h,guzai2,guzai3)-A(guzai1,guzai2,guzai3))/hdA_dguzai2=(A(guzai1,guzai2+h,guzai3)-A(guzai1,guzai2,guzai3))/hdA_dguzai3=(A(guzai1,guzai2,guzai3+h)-A(guzai1,guzai2,guzai3))/he_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=eph=0.001X=c(4,3,6)for(kin1:ncol(ep)){for(jin1:nrow(ep)){if(k==j){vector=X;vector_h=vectorvector_h[j]=vector_h[j]+hep[j,k]=(u(vector_h,k)-u(vector,k))/h}else{vector1=X;vector_h1=vector1vector_h1[j]=vector_h1[j]+hvector2=X;vector_h2=vector2vector_h2[k]=vector_h2[k]+hep[j,k]=((u(vector_h1,k)-u(vector1,k))/h+(u(vector_h2,j)-u(vector2,j))/h)/2abc[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%*%Xdu[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])