next up previous
Next: 5.3 二日目の晩 Up: 5 戦い Previous: 5.1 初日

5.2 二日目

これを書いているのは実は三日目なのだが、 大変な問題に関わってしまったと思い始めている。

初日ので大体話は済んだと思っていたのだが、 I君は別の計算を始めていた。 そもそもこの計算は Galerkin 法に基づいているのだが、 これまで使ってきた基底関数は境界条件を満足していないので、 あまりうまくないのだという。 境界条件を満たす基底関数を使うと収束が良くなる、と本に書いてあった。 ところがそこには基底関数が 2 つしか書かれていなかった:

$\displaystyle w_1(x)=x^2(6-4x+x^2), \quad w_2(x)=x^3(10-10x+3x^2) $

そこで基底関数の一般形を求めることが問題。 はっきりとは書いていないが、加藤先生の頭にあったのは

$\displaystyle w_n(x)=x^n\times( \hbox{$x$ の 2 次式で $1$ における 2 階、3 階の微分係数が $0$ になるもの}) $

ということらしい。こういうものを具体的に求めるのは高校数学の問題だが、 面倒だね。I君に紙の上で計算させたが、時間がかかりそうだったので、 Mathematica を持ち出して計算を始めてみた。 実際にどういう計算をしたのかは、 ログに残っているが、エッセンスをプログラムにまとめると、 
u=x^(n+1)(x^2+a x+b)
DDu=D[u,{x,2}]
u2=Simplify[DDu /. x->1]
DDDu=D[u,{x,3}]
u3=Simplify[DDDu /. x->1]
sol=Simplify[Solve[{u2==0,u3==0},{a,b}]]
ということである。その結果は 
oyabun% math
Mathematica 2.0 for SPARC
Copyright 1988-91 Wolfram Research, Inc.
 -- X11 windows graphics initialized --

In[1]:= << ishibashi1.m

                                           2
               -2 (3 + n)       6 + 5 n + n
Out[1]= {{a -> ----------, b -> ------------}}
                 1 + n                  2
                                   n + n

In[2]:=
となる。そこで

$\displaystyle w_n(x)=x^{n+1}[n(n+1)x^2-2n(n+3)x+(n+2)(n+3)] $

という基底関数の一般形が求まる。 これは $ n=1,2$ の時に加藤先生のあげたものと一致する ($ n=1$ の時は加藤先生の関数を $ 2$ 倍したもの)。

僕がこの計算を終える頃にはI君も大体計算を終えていて、 同じ結果を得た (まあ間違っていない限り同じになるのは当たり前だけど)。

これで話が済んだわけではない。 この基底関数を用いて行列 $ A$ , $ B$ を計算しないといけない。実は

$\displaystyle a_{ij}=(w''_i, w''_j), \quad b_{ij}=(w_i, w_j) $

ということだそうだ。ここで $ (,)$ は区間 $ (0,1)$ 上の内積を表す。

あまり手でやりたくない計算だ、 ということで Mathematica との格闘が始まった。 これも一部始終がログに残っているが、エッセンスをまとめると、 

%
u=x^(i+1)((i+1)i x^2-2i(i+3)x+(i+2)(i+3))
v=x^(j+1)((j+1)j x^2-2j(j+3)x+(j+2)(j+3))
intuv=Integrate[u v, {x,0,1}]
uv=Simplify[intuv]
Table[uv, {i,5}, {j,5}]
FortranForm[uv]
Print[Table[uv, {i,15}, {j,15}]]
Du=Simplify[D[u,{x,2}]]
Print["Du=",Du]
Dv=Simplify[D[v,{x,2}]]
Print["Dv=",Dv]
DuDv = Simplify[Du Dv]
Print["DuDv=",DuDv]
tmp=Integrate[DuDv /. i+j-2->m, {x,0,1}]
tmp=Simplify[tmp]
u2v2 = tmp/. m->i+j-2
Print["u2v2=", u2v2]
FortranForm[u2v2]
Print[Table[u2v2, {i,15}, {j,15}]]
ということになる。 

oyabun% math
Mathematica 2.0 for SPARC
Copyright 1988-91 Wolfram Research, Inc.
 -- X11 windows graphics initialized --

In[1]:= << foo.m
   1 + i                                               2
u=x      ((2 + i) (3 + i) - 2 i (3 + i) x + i (1 + i) x )
   1 + j                                               2
v=x      ((2 + j) (3 + j) - 2 j (3 + j) x + j (1 + j) x )
                              2        3      4
Bij=(8 (3780 + 3564 i + 1209 i  + 174 i  + 9 i  + 3564 j + 2505 i j +

              2         3           2          2       2  2        3
>        565 i  j + 40 i  j + 1209 j  + 565 i j  + 65 i  j  + 174 j  +

               3      4
>        40 i j  + 9 j )) /

>    ((3 + i + j) (4 + i + j) (5 + i + j) (6 + i + j) (7 + i + j))
TeXForm={{8\,\left( 3780 + 3564\,i + 1209\,{i^2} + 174\,{i^3} + 9\,{i^4} +

>         3564\,j + 2505\,i\,j + 565\,{i^2}\,j + 40\,{i^3}\,j + 1209\,{j^2} +

>         565\,i\,{j^2} + 65\,{i^2}\,{j^2} + 174\,{j^3} + 40\,i\,{j^3} +

>         9\,{j^4} \right) }\over

>     {\left( 3 + i + j \right) \,\left( 4 + i + j \right) \,

>       \left( 5 + i + j \right) \,\left( 6 + i + j \right) \,

>       \left( 7 + i + j \right) }}
FortranForm=8*(3780 + 3564*i + 1209*i**2 + 174*i**3 + 9*i**4 + 3564*j +

>       2505*i*j + 565*i**2*j + 40*i**3*j + 1209*j**2 + 565*i*j**2 +

>       65*i**2*j**2 + 174*j**3 + 40*i*j**3 + 9*j**4)/

>    ((3 + i + j)*(4 + i + j)*(5 + i + j)*(6 + i + j)*(7 + i + j))
                    2    3          2  -1 + i
Du=i (6 + 11 i + 6 i  + i ) (-1 + x)  x
                    2    3          2  -1 + j
Dv=j (6 + 11 j + 6 j  + j ) (-1 + x)  x
                      2    3                   2    3          4  -2 + i + j
DuDv=i (6 + 11 i + 6 i  + i ) j (6 + 11 j + 6 j  + j ) (-1 + x)  x
                         2    3                   2    3
Aij=(24 i (6 + 11 i + 6 i  + i ) j (6 + 11 j + 6 j  + j )) /

                                               2                  3
>    (120 + 274 (-2 + i + j) + 225 (-2 + i + j)  + 85 (-2 + i + j)  +

                      4               5
>      15 (-2 + i + j)  + (-2 + i + j) )
TeXForm={{24\,i\,\left( 6 + 11\,i + 6\,{i^2} + {i^3} \right) \,j\,

>       \left( 6 + 11\,j + 6\,{j^2} + {j^3} \right) }\over

>     {120 + 274\,\left( -2 + i + j \right)  +

>       225\,{{\left( -2 + i + j \right) }^2} +

>       85\,{{\left( -2 + i + j \right) }^3} +

>       15\,{{\left( -2 + i + j \right) }^4} +

>       {{\left( -2 + i + j \right) }^5}}}
FortranForm=24*i*(6 + 11*i + 6*i**2 + i**3)*j*(6 + 11*j + 6*j**2 + j**3)/

>    (120 + 274*(-2 + i + j) + 225*(-2 + i + j)**2 + 85*(-2 + i + j)**3 +

>      15*(-2 + i + j)**4 + (-2 + i + j)**5)
  416  2608  8636  1123  13384    2608  5216  227  904  17540
{{---, ----, ----, ----, -----}, {----, ----, ---, ---, -----},
  45   315   1155  165   2145     315   693   33   143  3003

     8636  227  908  5890  5480    1123  904  5890  2356  401
>   {----, ---, ---, ----, ----}, {----, ---, ----, ----, ---},
     1155  33   143  1001  1001    165   143  1001  429   78

     13384  17540  5480  401  3208
>   {-----, -----, ----, ---, ----}}
     2145   3003   1001  78   663
  576      576                960  1080
{{---, 96, ---, 72, 64}, {96, ---, ----, 160, 160},
   5        7                  7    7

     576  1080  1440       2880                  3360  3920
>   {---, ----, ----, 240, ----}, {72, 160, 240, ----, ----},
      7    7     7          11                    11    11

              2880  3920  62720
>   {64, 160, ----, ----, -----}}
               11    11    143

In[2]:= uv
後はI君に任せよう。Mathematica 万歳。

next up previous
Next: 5.3 二日目の晩 Up: 5 戦い Previous: 5.1 初日
桂田 祐史
2011-11-14