Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationEventAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
0
Help us understand the problem. What are the problem?

posted at

updated at

髪の毛の長さの時間変化

目的

髪の毛は切らなければ伸び続けるというのは、常識である。しかし、髪の毛が一定確率で抜けるのであれば、髪の毛の平均の長さは限界があるのではないだろうか。本紙では、数値シミュレーションを用いて髪の毛の長さの時間経過をシミュレーションし、その特徴を調べる。

方法

flowchart.png

髪の毛の本数は100000本1、ある髪が1日に髪の毛が抜ける確率を$\frac{1日に抜ける髪の本数}{髪の本数}$とし、1日に抜ける本数を100本2とした。また、髪の毛は抜けると同時に新しく生えるとし、髪全体の本数は不変であるとした。どの髪の毛も0.35 mm/day3で伸びるとする。髪がまったくない状態を初期条件にした。

結果

髪の毛の平均の長さは以下のようになった。

stat_ave.png

3000日前後から平均の髪の毛の長さが350 mmほどで安定してる。

髪の毛の長さのヒストグラムは以下のようになった。

hist.gif

多少の擾乱はあるが、指数分布に安定している。

考察

数値実験の結果から、髪の毛の平均の長さは無限に長くならないことが分かった。この結果を解析的に導けないか試してみる。

髪の毛の長さを$x$の本数を$h$とし、髪の毛の長さ分布の時間変化が成長と、抜けることによる影響しかないと仮定すると、
$$
h(x,t + \Delta t) \Delta x
= h(x,t) \Delta x
+ h\left(x - \frac{\Delta x}{2}- g \Delta t , t \right) g \Delta t
- h\left(x + \frac{\Delta x}{2}- g \Delta t , t \right) g \Delta t
- p \Delta t \cdot h(x,t) \Delta x
$$
と書ける。ただし、$g$は髪の毛の成長速度、$p$は髪の毛の抜ける確率だとする。$\Delta x \to 0, \Delta t \to 0$とすると、
$$
\frac{h(x, t + \Delta t) - h(x,t)}{\Delta t} =
- g \frac{ h\left(x + \frac{\Delta x}{2}- g \Delta t , t \right)

- h\left(x - \frac{\Delta x}{2}- g \Delta t , t \right) }{\Delta x}
- p h(x,t) \
\frac{\partial h}{\partial t}
= - g \frac{\partial h}{\partial x} - ph
$$
と微分方程式を導くことができる。ちなみにこの偏微分方程式の解は任意関数$h_0(x)$を用いて
$$
h(x,t) = h_0(x-gt)e^{-pt}
$$
と書けるなんで。定常$\frac{\partial h}{\partial t}=0$を仮定すると、
$$
g \frac{\partial h}{\partial x} = -ph \
\frac{1}{h} dh = - \frac{p}{g} dx \
\int \frac{1}{h} dh = - \frac{p}{g} \int dx \
\ln h = - \frac{p}{g} x + C \
h = A \exp \bigg( - \frac{p}{g} x \bigg)
$$
$h$の正規性$\int_0^\infty h(x) dx = 1$を用いると、
$$
\int_0^\infty A \exp{\left( - \frac{p}{g} x \right)} dx \
= A \left[ - \frac{g}{p} \exp{ \left( - \frac{p}{g} x \right)} \right]_0^\infty \
= A \frac{g}{p} = 1 \
A = \frac{p}{g}
$$
したがって、十分時間がたった髪の毛の長さの分布は、
$$
h(x) = \frac{p}{g} \exp{\left( - \frac{p}{g} x \right)}
$$
となる。ここで、$p = 10^{-3}$ [day^-1^]、$g = 0.35$ [mm day^-1^]を用いて、数値解に重ねると、

hist10000.png

となった。したがって、定常解は
$$
g \frac{\partial h}{\partial x} = -ph
$$
という微分方程式に従うと言える。

付録

fortranのコードを付けておきます。スレッド並列化したので、多少早く実行できます。

 program hair

   !$ use omp_lib

   implicit none

   integer, parameter :: time_max = 100*100
   integer, parameter :: hair_num = 100000 
   integer, parameter :: escape_num = 100  ! [/day]
   integer, parameter :: output_step = 100 
   real(8), parameter :: speed_hair = 0.35 ! [mm/day] 

   real(8), dimension(1:hair_num) :: length_hair

   real(8) :: prob_escape
   real(8) :: random
   real(8) :: total_length, total_length_2
   real(8) :: time_start, time_end 
   integer :: escaped_hair
   integer :: time_i, hair_i


   ! initialize

!   prob_escape = dble( escape_num ) / dble( hair_num )
   prob_escape = 1.d-3 
   length_hair(:) = 0.d0

   print *, "-----------------------------------"
   print *, " hair_num     = ", hair_num 
   print *, " speed_hair   = ", speed_hair
   print *, " prob_escape  = ", prob_escape
   print *, "-----------------------------------"

   call output( 0, hair_num, length_hair )

   open(11, file='./output/stat.txt', status='replace')

   !$ time_start = omp_get_wtime()
   do time_i = 1, time_max

     escaped_hair = 0  
     total_length   = 0.d0 
     total_length_2 = 0.d0 

     !$OMP parallel do private(random) reduction(+:escaped_hair) reduction(+:total_length) reduction(+:total_length_2)
     do hair_i = 1, hair_num

       call random_number( random )

       if ( random < prob_escape ) then
       ! escape hair

         length_hair(hair_i) = 0.d0      
         escaped_hair = escaped_hair + 1

       end if

       length_hair(hair_i) = length_hair(hair_i) + speed_hair
       total_length   = total_length   + length_hair(hair_i)   
       total_length_2 = total_length_2 + length_hair(hair_i)**2   

     end do
     !$OMP end parallel do

     if ( mod( time_i, output_step ) == 0 ) then

       print *, time_i, total_length/dble(hair_num)

       call output( time_i, hair_num, length_hair )

     end if

     write(11, '(I5.5, 3X, I3.3, 3X, 3(f15.5, 3X))') time_i, &
                 escaped_hair,                     &
                 total_length/dble(hair_num),      &
                 total_length_2/dble(hair_num)-(total_length/dble(hair_num))**2, &
                 maxval(length_hair(:))         

   end do  
   !$ time_end = omp_get_wtime()
   !$ print *, "time = ", time_end - time_start

   close(11)


   contains

   subroutine output( time_i, hair_num, length_hair )

     implicit none

     integer,                        intent(in) :: time_i
     integer,                        intent(in) :: hair_num
     real(8), dimension(1:hair_num), intent(in) :: length_hair

     integer :: hair_i
     character(128) :: filename

     write(filename, '("./output/length/length"I5.5".txt")') time_i

     open(10, file=filename, status='replace' )

       do hair_i = 1, hair_num

         write(10,*) hair_i, length_hair(hair_i)

       end do

     close(10)

   end subroutine output

 end program hair

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
0
Help us understand the problem. What are the problem?