DE積分 振動型積分 cos版

\int_0^\infty f(x)\cos\omg x dx
f(x)=-exp(x)

module exp_lib
implicit none
real(8) pi

contains

subroutine calf(f,x)
real(8),intent(in) :: x
real(8),intent(out) :: f
f=exp(-x)
end subroutine calf

subroutine de2_int(sum,em,omg,n,calf)
integer,intent(in) :: n
real(8),intent(in) :: em,omg
real(8),intent(out) :: sum
real(8) h,t,e6,phi,dphi,x,f
integer i
h=pi/(em*omg)
t=0.d0+h*0.5d0
!!$ e6=1.d0
!!$ phi=1.d0/6.d0
!!$ dphi=1.d0/2.d0
e6=exp(-6.d0*sinh(t))
phi=t/(1.d0-e6)
dphi=(1.d0+(-6.d0*t*cosh(t)-1.d0)*e6)/(1.d0-e6)**2
x=em*phi
call calf(f,x)
sum=sum+f*cos(omg*x)*dphi
do i=1,n
t=i*h+h*0.5d0
e6=exp(-6.d0*sinh(t))
phi=t/(1.d0-e6)
dphi=(1.d0+(-6.d0*t*cosh(t)-1.d0)*e6)/(1.d0-e6)**2
x=em*phi
call calf(f,x)
sum=sum+f*cos(omg*x)*dphi
t=-t+h*0.5d0
e6=exp(-6.d0*sinh(t))
phi=t/(1.d0-e6)
dphi=(1.d0+(-6.d0*t*cosh(t)-1.d0)*e6)/(1.d0-e6)**2
x=em*phi
call calf(f,x)
sum=sum+f*cos(omg*x)*dphi
end do
sum=sum*em*h
end subroutine de2_int

end module exp_lib

program exp
use exp_lib
implicit none
external calf
integer n
real(8) em,omg,sum
pi=acos(-1.d0)
sum=0.d0
em=100.d0
omg=1.d0
n=75
call de2_int(sum,em,omg,n,calf)
write(*,*) sum
end program exp

DE積分 オイラーの定数

Fortranを使って振動型の積分
\int_0^\infty f(x)\sim\omg x dx
hが刻み幅
nは打ち切る項数
omgは周期
Mはhとomgから決まる数
詳しくは数値解析、森正武、共立出版p214


module euler_lib
implicit none
real(8) pi

contains

subroutine calf(f,x)
real(8),intent(in) :: x
real(8),intent(out) :: f
if(x == 0.d0) then
f=0.d0
else
f=log(x)
end if
end subroutine calf

subroutine de_int(sum,em,omg,n,calf)
integer,intent(in) :: n
real(8),intent(in) :: em,omg
real(8),intent(out) :: sum
real(8) h,t,e6,phi,dphi,x,f
integer i
h=pi/(em*omg)
t=0.d0
e6=1.d0
phi=1.d0/6.d0
dphi=1.d0/2.d0
x=em*phi
call calf(f,x)
sum=sum+f*sin(omg*x)*dphi
do i=1,n
t=i*h
e6=exp(-6.d0*sinh(t))
phi=t/(1.d0-e6)
dphi=(1.d0+(-6.d0*t*cosh(t)-1.d0)*e6)/(1.d0-e6)**2
x=em*phi
call calf(f,x)
sum=sum+f*sin(omg*x)*dphi
t=-t
e6=exp(-6.d0*sinh(t))
phi=t/(1.d0-e6)
dphi=(1.d0+(-6.d0*t*cosh(t)-1.d0)*e6)/(1.d0-e6)**2
x=em*phi
call calf(f,x)
sum=sum+f*sin(omg*x)*dphi
end do
sum=sum*em*h
end subroutine de_int

end module euler_lib

program euler
use euler_lib
implicit none
external calf
integer n
real(8) em,omg,sum
pi=acos(-1.d0)
sum=0.d0
em=50.d0
omg=1.d0
n=75
call de_int(sum,em,omg,n,calf)
write(*,*) sum
end program euler

円孔1/4対称、四角形メッシュ

トンネルのようなモデル
赤はGUIで操作した部分
emacsで編集、GUIで編集を繰り返して作成

l = 2 ;

Point(10) = {0, 0, 0, l};
Point(20) = {1, 0, 0, l} ;
Point(30) = {0, 1, 0, l} ;
Point(40) = {2, 0, 0, l} ;
Point(50) = {0, 2, 0, l} ;
Point(60) = {Sin(Pi/4),Sin(Pi/4),0,l} ;
Point(70) = {2*Sin(Pi/4),2*Sin(Pi/4),0,l} ;
Point(100)={3, 0, 0, l} ;
Point(110)={3, 3, 0, l} ;
Point(120)={0, 3, 0, l} ;


Circle(1) = {20,10,60};
Circle(2) = {60,10,30};
Circle(3) = {40,10,70};
Circle(4) = {70,10,50};
Line(5) = {20,40};
Line(6) = {40,100};
Line(7) = {100,110};
Line(8) = {110,120};
Line(9) = {120,50};
Line(10) = {50,30};
Line(11) = {60,70};
Line(12) = {70,110};

Transfinite Line{1,3,2,4,8,7}=5+1;
Transfinite Line{5,11,10}=2+1;
Transfinite Line{6,12,9}=4+1;

Line Loop(13) = {1,11,-3,-5};
Plane Surface(14) = {13};
Transfinite Surface{14}={20,40,70,60};
Recombine Surface{14};

Line Loop(15) = {2,-10,-4,-11};
Plane Surface(16) = {15};
Transfinite Surface{16}={60,70,50,30};
Recombine Surface{16};

Line Loop(17) = {12,8,9,-4};
Plane Surface(18) = {17};
Transfinite Surface{18}={70,50,120,110};
Recombine Surface{18};

Line Loop(19) = {7,-12,-3,6};
Plane Surface(20) = {19};
Transfinite Surface{20}={40,100,110,70};
Recombine Surface{20};

Physical Surface(21) = {20,14,16,18};

バームクーヘンを四角形要素で

l = 2 ;

Point(1) = {0, 0, 0, l};
Point(2) = {1, 0, 0, l} ;
Point(3) = {0, 1, 0, l} ;
Point(4) = {2, 0, 0, l} ;
Point(5) = {0, 2, 0, l} ;

Circle(1001)={2,1,3};
Circle(1002)={5,1,4};
Line(1003)={4,2};
Line(1004)={3,5};

Line Loop(1101)={1001,1004,1002,1003};
Plane Surface(2001)={1101};

Transfinite Line{1001,1002}=11;
Transfinite Line{1003,1004}=3;
Transfinite Surface{2001}={2,3,5,4};
Recombine Surface {2001};

MySurface = 100;
Physical Surface(MySurface) = {3001} ;

正方形を四角形要素で分割

l = 2 ;

Point(1) = {0, 0, 0, l};
Point(2) = {1, 0, 0, l} ;
Point(3) = {1, 1, 0, l} ;
Point(4) = {0, 1, 0, l} ;

Line(11) = {1,2} ;
Line(12) = {2,3} ;
Line(13) = {3,4} ;
Line(14) = {4,1} ;

Line Loop(5) = {11,12,13,14} ;

Plane Surface(6) = {-5} ;
Transfinite Line{11,12,13,14}=10;
Transfinite Surface{6}={1,2,3,4};
Recombine Surface {6};

// Consequently, two punctual elements will be saved in the output
// files, both with the region number 1. The mechanism is identical
// for line or surface elements:

// Physical Line(10) = {1,2,4} ;

MySurface = 100;
Physical Surface(MySurface) = {6} ;