« 2025年6月 | トップページ

2025年7月

2025年7月11日 (金)

垂水の皇大神社、猿田彦大神、倉稲魂大神、白龍大神、黒龍大神が一か所に集まっているところでお参り。

垂水駅から細い道を行くと現れます。小さな社が3つあります。

20241123-160156 20241123-160213 20241123-160216 20241123-160222

2025年7月10日 (木)

Interface2025年8月号Pythonで体験!はじめての暗号を買った。上杉暗号からRSA、AES、DHなど、特に楕円曲線暗号についてはコードも実際に動かすところまで詳しくかかれていた。耐量子暗号や聞いたことなかったY-00暗号や関数型暗号も記載。

暗号は昔、結城浩さんの暗号技術入門を読んだが、そこから全然アップデートしてないので読んでみた。

こういうのも作ったし。

RSA暗号の検証(NHK笑わない数学、暗号理論より)ができる自作式をカシオの高精度計算サイトkeisan.casio.jpにUP!暗号化と復号化ができます。アルゴリズムは結城浩さんの暗号技術入門を参考にしました。冪剰余、最小公倍数、拡張ユークリッドの互除法も詰め込んだ。

20250709-174426

アマゾンリンク:https://amzn.to/4nBlzpw

内容はとても多岐にわたっていて、セキュリティの考え方から上杉暗号、シーザー暗号、たぬき暗号!などから始まり(エニグマを模したコードまで)、RSA、AES、そして楕円曲線暗号についてはかなりの紙面を割いて最初の考え方からコード、実際にラズパイpico 2 wで動かすところまで乗っていて面白かった。

こんなにコンピュータの雑誌で楕円曲線が書かれているのを見たのは初めて。

それ以外にも形式検証ツールや耐量子暗号、関数型暗号、Y-00暗号と話題が豊富。TLSの実装まで。

これでだいぶ暗号についてはアップデートした気がする。

 

2025年7月 9日 (水)

天田神社(あまたじんじゃ、大阪の交野市)でお参り。

ここは初めて来た。京阪の河内森駅からすぐです。

20241221-150228

20241221-150020

20241221-150029

20241221-150030

20241221-150036

 

高周波・RFニュース 2025年7月8日 NordicとSercommのセルラーIoTモジュール、iFixitがFairphone 6を分解、スコアは10/10、RCR wireless newsのウェビナー2件(6GとIndustry4.0)、SEMCOが高耐圧C0G MLCCを車載急速充電に提案

・NordicとSercommのセルラーIoTモジュール

Nordic Semiconductor and Sercomm Corporation partner on cellular IoT module designed for applications demanding reliable, power-saving performance

202507091

・iFixitがFairphone 6を分解、スコアは10/10

ポゴピン接続が多いのはなるほど。

Fairphone 6 Teardown: Proof Phones Don’t Have to Be Disposable

202507092

・RCR wireless newsのウェビナー2件(6GとIndustry4.0)

https://content.rcrwireless.com/webinar-6g-standards-spectrum-use-cases

https://content.rcrwireless.com/powering-logistics-4.0-webinar

202507093

・SEMCOが高耐圧C0G MLCCを車載急速充電に提案

Proposal for High-Voltage MLCCs in Alignment with xEV Fast-Charging Trends

202507094

2025年7月 8日 (火)

Google ColabのJulia言語で1次元のGray-Scottモデル(∂u/∂t=u²v-(F+k)u+Du∂²u/∂x²,∂v/∂t=-u²v+F(1-v)+Dv∂²v/∂x²)を計算してパルスが次々分裂する様子を見る。空間6次の差分、時間8次のルンゲクッタ法で計算。

今回はGray-Scott方程式を計算。2次元だと動物の模様のようなものができて面白いが今回は1次元でやってみよう。

∂u/∂t=u²v-(F+k)u+Du∂²u/∂x²
∂v/∂t=-u²v+F(1-v)+Dv∂²v/∂x²

コードはこんな感じ。

using Plots
using Printf

function main()
    n = 200 #空間分割数
    m = 2000000 #時間分割数

    #u, vの初期化
    u = zeros(Float64, n + 6)
    v = zeros(Float64, n + 6)
    ud = zeros(Float64, n + 6)
    vd = zeros(Float64, n + 6)


    #物理パラメータ
    xmax = 1.6
    xmin = 0.0
    t = 0.0

    du = 0.00001
    dv = 0.00002
    Fc = 0.04
    kc = 0.06075

    dt = 0.01
    dx = (xmax - xmin) / (n - 3)
    x = zeros(Float64, n + 6)

    #表示用
    results = []
    ndiv = 2000

    #8次のルンゲクッタ法の係数
    a = zeros(Float64, 13)
    b = zeros(Float64, 13, 13)
    c = zeros(Float64, 13)
    ku = zeros(Float64, n, 13)
    kv = zeros(Float64, n, 13)
    a[1] = 0.0
    a[2] = 2.0 / 27.0
    a[3] = 1.0 / 9.0
    a[4] = 1.0 / 6.0
    a[5] = 5.0 / 12.0
    a[6] = 0.5
    a[7] = 5.0 / 6.0
    a[8] = 1.0 / 6.0
    a[9] = 2.0 / 3.0
    a[10] = 1.0 / 3.0
    a[11] = 1.0
    a[12] = 0.0
    a[13] = 1

    b[2, 1] = 2.0 / 27.0

    b[3, 1] = 1.0 / 36
    b[3, 2] = 1.0 / 12.0

    b[4, 1] = 1.0 / 24.0
    b[4, 3] = 1.0 / 8.0

    b[5, 1] = 5.0 / 12.0
    b[5, 3] = -25.0 / 16.0
    b[5, 4] = 25.0 / 16.0

    b[6, 1] = 1.0 / 20.0
    b[6, 4] = 1.0 / 4.0
    b[6, 5] = 1.0 / 5.0

    b[7, 1] = -25.0 / 108.0
    b[7, 4] = 125.0 / 108.0
    b[7, 5] = -65.0 / 27.0
    b[7, 6] = 125.0 / 54.0

    b[8, 1] = 31.0 / 300.0
    b[8, 5] = 61.0 / 225.0
    b[8, 6] = -2.0 / 9.0
    b[8, 7] = 13.0 / 900.0

    b[9, 1] = 2.0
    b[9, 4] = -53.0 / 6.0
    b[9, 5] = 704.0 / 45.0
    b[9, 6] = -107.0 / 9.0
    b[9, 7] = 67.0 / 90.0
    b[9, 8] = 3.0

    b[10, 1] = -91.0 / 108.0
    b[10, 4] = 23.0 / 108.0
    b[10, 5] = -976.0 / 135.0
    b[10, 6] = 311.0 / 54.0
    b[10, 7] = -19.0 / 60.0
    b[10, 8] = 17.0 / 6.0
    b[10, 9] = -1.0 / 12.0

    b[11, 1] = 2383.0 / 4100.0
    b[11, 4] = -341.0 / 164.0
    b[11, 5] = 4496.0 / 1025.0
    b[11, 6] = -301.0 / 82.0
    b[11, 7] = 2133.0 / 4100.0
    b[11, 8] = 45.0 / 82.0
    b[11, 9] = 45.0 / 164.0
    b[11, 10] = 18.0 / 41.0

    b[12, 1] = 3.0 / 205.0
    b[12, 6] = -6.0 / 41.0
    b[12, 7] = -3.0 / 205.0
    b[12, 8] = -3.0 / 41.0
    b[12, 9] = 3.0 / 41.0
    b[12, 10] = 6.0 / 41.0

    b[13, 1] = -1777.0 / 4100.0
    b[13, 4] = -341.0 / 164.0
    b[13, 5] = 4496.0 / 1025.0
    b[13, 6] = -289.0 / 82.0
    b[13, 7] = 2193.0 / 4100.0
    b[13, 8] = 51.0 / 82.0
    b[13, 9] = 33.0 / 164.0
    b[13, 10] = 12.0 / 41.0
    b[13, 12] = 1.0

    c[6] = 34.0 / 105.0
    c[7] = 9.0 / 35.0
    c[8] = 9.0 / 35.0
    c[9] = 9.0 / 280.0
    c[10] = 9.0 / 280.0
    c[12] = 41.0 / 840.0
    c[13] = 41.0 / 840.0


    #初期値
    for i in 1:(n + 6)
        x[i] = xmin + (i - 3) * dx
    end
    for i in 4:(n + 3)
        u[i] = exp(-(x[i] - xmax / 2.0) ^ 2 * 1000.0)
        v[i] = 1.0 - u[i]
    end

    for i in 1:m

        # 8次のルンゲクッタ法計算
         @inbounds for ii in 4:(n + 3)
            ku[ii - 3, 1] = f(ii, t + dt, u, v, du, dv, Fc, kc, dx, n) * dt
            kv[ii - 3, 1] = g(ii, t + dt, u, v, du, dv, Fc, kc, dx, n) * dt
        end

        @inbounds for j in 2:13

            @inbounds @simd for ii in 4:(n + 3)
                ud[ii] = u[ii]
                vd[ii] = v[ii]
            end

            @inbounds for k in 1:(j - 1)
                @inbounds @simd for ii in 2:(n + 1)
                    ud[ii] = ud[ii] + b[j, k] * ku[ii - 3, k]
                    vd[ii] = vd[ii] + b[j, k] * kv[ii - 3, k]
                end
            end

            @inbounds for ii in 4:(n + 3)
                ku[ii - 3, j] = f(ii, t + a[j] * dt, ud, vd, du, dv, Fc, kc, dx, n) * dt
                kv[ii - 3, j] = g(ii, t + a[j] * dt, ud, vd, du, dv, Fc, kc, dx, n) * dt
            end
        end

        @inbounds @simd for j in 1:13
            for ii in 2:(n + 1)
                u[ii] = u[ii] + c[j] * ku[ii - 3, j]
                v[ii] = v[ii] + c[j] * kv[ii - 3, j]
            end
        end

        t = t + dt

        # ndivステップごとに結果(u)を配列に格納
        if i % ndiv == 0
            push!(results, copy(u))
        end

    end
    # 計算結果をアニメーションで表示
    anim = @animate for i in 1:length(results)
        plot(x, results[i], ylim = (0.0, 0.5), linewidth = 4,title="Gray-Scott model", label="", xlabel="x", ylabel="u", size=(900,500))
    end
    gif(anim, "GrayScott.gif", fps = 30)
end

function f(i, t, u, v, du, dv, Fc, kc, dx, n)
    u[1] = u[n + 1]
    u[2] = u[n + 2]
    u[3] = u[n + 3]
    u[n + 4] = u[4]
    u[n + 5] = u[5]
    u[n + 6] = u[6]

    d2 = (2 * u[i + 3] - 27.0 * u[i + 2] + 270.0 * u[i + 1] - 490.0 * u[i] + 270.0 * u[i - 1] - 27.0 * u[i - 2] + 2 * u[i - 3]) / (180.0 * dx * dx)

    return du * d2 + u[i] * u[i] * v[i] - (Fc + kc) * u[i]
end

function g(i, t, u, v, du, dv, Fc, kc, dx, n)
    v[1] = v[n + 1]
    v[2] = v[n + 2]
    v[3] = v[n + 3]
    v[n + 4] = v[4]
    v[n + 5] = v[5]
    v[n + 6] = v[6]

    d2 = (2 * v[i + 3] - 27.0 * v[i + 2] + 270.0 * v[i + 1] - 490.0 * v[i] + 270.0 * v[i - 1] - 27.0 * v[i - 2] + 2 * v[i - 3]) / (180.0 * dx * dx)

    return dv * d2 - u[i] * u[i] * v[i] + Fc * (1.0 - v[i])
end

実行結果。パルスが次々分裂していくのが見える。

 

 

 

2025年7月 7日 (月)

廣田神社(西宮)でお参り。

西宮震災記念碑公園の後はこちらの廣田神社へ。かなり立派な神社でした。ただここにたどり着くのに逆方向の山側から入ったので疲れた…

20250521-135610

20250521-135351

20250521-135452

20250521-135505

20250521-135234

20250521-135334

2025年7月 6日 (日)

元祖豚丼屋TONTON(上新庄)で豚バラ丼大盛をいただく。思っていた以上に厚いバラ肉がたくさんで、甘めのタレもかかって満足度高し。

ずっと前から行ってみたかったTONTONにようやく行けた。

特盛にしようかと思ったが様子見でまずは大盛で。

20250504-123645

20250504-123650

大盛でも想像以上の枚数のバラ肉で、厚みも十分。脂身も適度に入っていてなかなか美味しかった。

タレは甘めで、後半ちょっと空きかけるが卓上に山椒と七味唐辛子があるので味変もOK。

これはまた行きたい。

 

2025年7月 5日 (土)

正倉院 THE SHOW@大阪歴史博物館へ行ってきた。蘭奢待の香りを再現して嗅がせてくれるもの、確かに上品でとてもいい香りだった。奈良の正倉院展には開催されるたびに行っているがものすごい混雑なので模造品でも近くでゆっくり観られるのはいい。

久しぶりに大阪歴史博物館へ。NHKの隣にある。

20250704-112937

20250704-115405

模造品とはいえよくできた正倉院の宝物が間近で観られるのはいい。

20250704-113447

20250704-113609

20250704-113629

そして蘭奢待。

20250704-114544

香りが再現されているということでぜひ嗅ぎたいということで来たのだった。

20250704-114555

並ぶのかな…と思っていたら20個の小型ガラス容器で体験できるようになっていて賢いなと思った。

で嗅いだ感想は、やっぱり上品で強いいい香りがした。

20250704-114602

2025年7月 4日 (金)

泉大津のシーパスパークをぶらぶら歩く。

ものすごくいい天気でした。

20250505-135620

20250505-135603

20250505-135535

20250505-135522 20250505-132346

この時はフラダンスをやっていた。

20250505-132232

20250505-132111

20250505-132102

 

 

2025年7月 3日 (木)

夙川公園で水車を見た。

残念ながら復元されたものではなかった…

20250521-131958

20250521-132003

2025年7月 2日 (水)

ローズマリーのあまき香り(島田荘司さん)を読んだ。20年前のニューヨークでバレエ上演中に伝説のバレリーナが幕合に密室で殺されるが、観客は最後まで踊りをみたという。20年後に御手洗潔が真相を解き明かす。ユダヤ、ナチス、ウクライナ、日本、など興味ある話題も豊富。

本当に久しぶりに島田荘司さんの御手洗潔シリーズを読んだ。

20250701-171505

アマゾンリンク:https://amzn.to/4ev8t9o

あらすじは
「世界中で人気を博す、生きる伝説のバレリーナ・クレスパンが密室で殺された。
1977年10月、ニューヨークのバレエシアターで上演された「スカボロゥの祭り」で主役を務めたクレスパン。
警察の調べによると、彼女は2幕と3幕の間の休憩時間の最中に、専用の控室で撲殺されたという。
しかし3幕以降も舞台は続行された。
さらに観客たちは、最後までクレスパンの踊りを見ていた、と言っていてーー?」

というもの。もう最初からわけがわからない展開。

しかし20年後の1997年、いつの間にかスウェーデンの大学で脳を研究している御手洗潔(ネジ式ザゼツキーでそういう描写があったとか)が、その事件を調べている友人のジャーナリストに頼まれて推理をしていく。スウェーデン時代なので石岡さんは残念ながらでてきません。

そして3つの全く関係ない、アメリカで起きた間抜けな事件を聞いた御手洗さんは真相に迫っていく…

トリックの一つはある別の作家さんのものすごく有名な作品(映像化されるそうだ)を思い出し、もう一つは島田さんの別の作品を思い出したが、これはトリック云々よりそのクレスパンの人生や、それにかかわるユダヤ、ナチスなどなどの背景が重要で、また今の戦争についての批判もあり島田さんらしい。

長い作品だが面白くて一気に読めた。

 

 

高周波・RFニュース 2025年7月2日 5G Americasが6Gに向けセンシングと通信ホワイトペーパー発行、KYOCERA AVXが3dBハイブリッドカプラ発表、TDKが車載薄膜インダクタ発表、Nordicが1次電池向けPMIC発表、ローデ・シュワルツの6GとAI/ML解説記事

・5G Americasが6Gに向けセンシングと通信ホワイトペーパー発行

5G Americas Highlights the Strategic Role of Integrated Sensing and Communication for 6G

202507021

・KYOCERA AVXが3dBハイブリッドカプラ発表

KYOCERA AVX Releases New 3dB Hybrid Couplers

202507022

・TDKが車載薄膜インダクタ発表

インダクタ: 車載電源回路用薄膜インダクタの開発と量産

202507023

・Nordicが1次電池向けPMIC発表

Nordic Semiconductor’s nPM2100 Power Management IC for primary battery applications enters mass production

202507024

・ローデ・シュワルツの6GとAI/ML解説記事

AI and ML for 6G networks

・DOCSIS4.0で14Gbps達成

CableLabs sees 14-Gig speeds at latest DOCSIS 4.0 interop

2025年7月 1日 (火)

Google ColabのJulia言語で時空カオスを生む蔵本・シバシンスキー方程式(Kuramoto–Sivashinsky equation)∂u/∂t+u∂u/∂x+∂²u/∂²+∂⁴/∂x⁴=0を空間6次の差分、時間は8次のルンゲクッタ法で計算して動画にしてみる。

今回は蔵本・シバシンスキー方程式。時間・空間的にカオスを生じる最も簡単な偏微分方程式と言われている。金平糖の方程式と呼ばれているとかいないとか。

∂u/∂t+u∂u/∂x+∂²u/∂²+∂⁴/∂x⁴=0

これは何回もExcel VBAで計算しているのでそれを移植した。

空間方向は6次の差分(4階微分が出てくるので他の微分もそれに合わせた)でこんな感じで書ける。

1階: (-u(x-3dx)/9 + u(x-2dx)-5u(x-dx)+5u(x+dx) -u(x+2dx)+u(x+3dx)/9)/(20dx/3)
2階: (2u(x+3dx)-27u(x+2dx)+270u(x+dx)-490u(x)+270u(x-dx)-27u(x-2dx)+2u(x-3dx))/(180dx^2)
3階: (-u(x+3dx)+8u(x+2dx)-13u(x+dx)+13u(x-dx)-8u(x-2dx)+u(x-3dx))/(8dx^3)
4階: (-u(x+3dx)+12u(x+2dx)-39u(x+dx)+56u(x)-39u(x-dx)+12u(x-2dx)-u(x-3dx))/(6dx^4)

時間方向はRunge Kutta Felhbergの8次を使う。

コードはこんな感じ。


using Plots
using Random
using Printf

function main()
    n = 200 #空間分割数
    m = 400000 #時間分割数

    u = zeros(Float64, n + 6)
    ud = zeros(Float64, n + 6)


    xmax = 50.0
    xmin = 0.0
    x = zeros(Float64, n + 6)
    D = 0.15

    t = 0.0
    dt = 0.0005
    dx = (xmax - xmin) / (n - 3)
    MT = MersenneTwister()
    Random.seed!(MT, 42)

    #表示用
    results = []
    ndiv = 500

    #8次のルンゲクッタ法の係数
    a = zeros(Float64, 13)
    b = zeros(Float64, 13, 13)
    c = zeros(Float64, 13)
    ku = zeros(Float64, n, 13)

    a[1] = 0.0
    a[2] = 2.0 / 27.0
    a[3] = 1.0 / 9.0
    a[4] = 1.0 / 6.0
    a[5] = 5.0 / 12.0
    a[6] = 0.5
    a[7] = 5.0 / 6.0
    a[8] = 1.0 / 6.0
    a[9] = 2.0 / 3.0
    a[10] = 1.0 / 3.0
    a[11] = 1.0
    a[12] = 0.0
    a[13] = 1

    b[2, 1] = 2.0 / 27.0

    b[3, 1] = 1.0 / 36
    b[3, 2] = 1.0 / 12.0

    b[4, 1] = 1.0 / 24.0
    b[4, 3] = 1.0 / 8.0

    b[5, 1] = 5.0 / 12.0
    b[5, 3] = -25.0 / 16.0
    b[5, 4] = 25.0 / 16.0

    b[6, 1] = 1.0 / 20.0
    b[6, 4] = 1.0 / 4.0
    b[6, 5] = 1.0 / 5.0

    b[7, 1] = -25.0 / 108.0
    b[7, 4] = 125.0 / 108.0
    b[7, 5] = -65.0 / 27.0
    b[7, 6] = 125.0 / 54.0

    b[8, 1] = 31.0 / 300.0
    b[8, 5] = 61.0 / 225.0
    b[8, 6] = -2.0 / 9.0
    b[8, 7] = 13.0 / 900.0

    b[9, 1] = 2.0
    b[9, 4] = -53.0 / 6.0
    b[9, 5] = 704.0 / 45.0
    b[9, 6] = -107.0 / 9.0
    b[9, 7] = 67.0 / 90.0
    b[9, 8] = 3.0

    b[10, 1] = -91.0 / 108.0
    b[10, 4] = 23.0 / 108.0
    b[10, 5] = -976.0 / 135.0
    b[10, 6] = 311.0 / 54.0
    b[10, 7] = -19.0 / 60.0
    b[10, 8] = 17.0 / 6.0
    b[10, 9] = -1.0 / 12.0

    b[11, 1] = 2383.0 / 4100.0
    b[11, 4] = -341.0 / 164.0
    b[11, 5] = 4496.0 / 1025.0
    b[11, 6] = -301.0 / 82.0
    b[11, 7] = 2133.0 / 4100.0
    b[11, 8] = 45.0 / 82.0
    b[11, 9] = 45.0 / 164.0
    b[11, 10] = 18.0 / 41.0

    b[12, 1] = 3.0 / 205.0
    b[12, 6] = -6.0 / 41.0
    b[12, 7] = -3.0 / 205.0
    b[12, 8] = -3.0 / 41.0
    b[12, 9] = 3.0 / 41.0
    b[12, 10] = 6.0 / 41.0

    b[13, 1] = -1777.0 / 4100.0
    b[13, 4] = -341.0 / 164.0
    b[13, 5] = 4496.0 / 1025.0
    b[13, 6] = -289.0 / 82.0
    b[13, 7] = 2193.0 / 4100.0
    b[13, 8] = 51.0 / 82.0
    b[13, 9] = 33.0 / 164.0
    b[13, 10] = 12.0 / 41.0
    b[13, 12] = 1.0

    c[6] = 34.0 / 105.0
    c[7] = 9.0 / 35.0
    c[8] = 9.0 / 35.0
    c[9] = 9.0 / 280.0
    c[10] = 9.0 / 280.0
    c[12] = 41.0 / 840.0
    c[13] = 41.0 / 840.0


    #初期値
    for i in 1:(n + 6)
        x[i] = xmin + (i - 3) * dx
    end
    for i in 4:(n + 3)
        u[i] = 0.1 * (2.0 * rand(MT) - 1.0)
    end
    u[1] = u[n + 1]
    u[2] = u[n + 2]
    u[3] = u[n + 3]
    u[n + 4] = u[4]
    u[n + 5] = u[5]
    u[n + 6] = u[6]

    for i in 1:m

        # 8次のルンゲクッタ法計算
         @inbounds for ii in 4:(n + 3)
            ku[ii - 3, 1] = fu(ii, t + dt, u, dx, n) * dt
          end

        @inbounds for j in 2:13

            @inbounds @simd for ii in 4:(n + 3)
                ud[ii] = u[ii]
            end

            @inbounds for k in 1:(j - 1)
                @inbounds @simd for ii in 4:(n + 3)
                    ud[ii] = ud[ii] + b[j, k] * ku[ii - 3, k]
                end
            end

            @inbounds for ii in 4:(n + 3)
                ku[ii - 3, j] = fu(ii, t + a[j] * dt, ud, dx, n) * dt

            end
        end

        @inbounds @simd for j in 1:13
            for ii in 4:(n + 3)
                u[ii] = u[ii] + c[j] * ku[ii - 3, j]
            end
        end
        t = t + dt

        # ndivステップごとに結果を配列に格納
        if i % ndiv == 0
            push!(results, copy(u))
        end

    end
    # 計算結果をアニメーションで表示
    anim = @animate for i in 1:length(results)
        plot(x, results[i], ylim = (-3, 3), linewidth = 3,title="Kuramoto Sivashinsky equation", label="", xlabel="x", ylabel="u", size=(900,500))
    end
    gif(anim, "Kuramoto_Sivashinsky.gif", fps = 30)
end

 function fu(i, t, u, dx, n)

    u[1] = u[n + 1]
    u[2] = u[n + 2]
    u[3] = u[n + 3]
    u[n + 4] = u[4]
    u[n + 5] = u[5]
    u[n + 6] = u[6]

    d1 = (u[i + 3] / 9.0 - u[i + 2] + 5.0 * u[i + 1] - 5.0 * u[i - 1] + u[i - 2] - u[i - 3] / 9.0) / (20.0 * dx / 3.0)
    d2 = (2.0 * u[i + 3] - 27.0 * u[i + 2] + 270.0 * u[i + 1] - 490.0 * u[i] + 270.0 * u[i - 1] - 27.0 * u[i - 2] + 2 * u[i - 3]) / (180.0 * dx * dx)
    d4 = (-u[i + 3] + 12.0 * u[i + 2] - 39.0 * u[i + 1] + 56.0 * u[i] - 39.0 * u[i - 1] + 12.0 * u[i - 2] - u[i - 3]) / (6.0 * dx * dx * dx * dx)

    return -u[i] * d1 - d2 - d4

 end

main()で実行させるとこうなる。

全く予想できない動きをしている。

« 2025年6月 | トップページ

最近の記事

2025年7月
    1 2 3 4 5
6 7 8 9 10 11 12
13 14 15 16 17 18 19
20 21 22 23 24 25 26
27 28 29 30 31    

最近のコメント

無料ブログはココログ
フォト