

機械系の学生ですが、課題の一部にプログラミングがあり以下のコードでポアソン方程式を解きCSVファイルへ計算結果を出力する方法を模索していましたが数日進まないので質問させていただきます。プログラミング自体初心者なのでなるべく初心者でもわかるように説明していただければ幸いです。ご容赦ください。Google
Colaboratoryを使用しています。  
質問内容  
エラーは出ませんが実行すると一列の数字が出力され３つの列になりません。そもそもx単体で出力すると何も出てこなかったです。グラフを見る限り計算自体は合っていると思うのでどうにかできないかと思いここに来ました。[![計算を一回で出しているので思った通りの形になってます](https://i.stack.imgur.com/1DRHQ.png)](https://i.stack.imgur.com/1DRHQ.png)  
[![実行結果です](https://i.stack.imgur.com/sLMm0.png)](https://i.stack.imgur.com/sLMm0.png)  
一番いいのはエクセルにそのまま出力できるのが理想ですが難しそうなのでとりあえず出力できることだけを目標に行ってます。

[code]

    from logging import fatal
    from matplotlib import pyplot as plt
    from matplotlib import cm
    from mpl_toolkits.mplot3d import Axes3D
    import numpy as np
    import itertools
    import csv
    
    # 条件設定
    nx = 100
    ny = 100
    nt  = 1
    xmin = 0
    xmax = 2
    ymin = 0
    ymax = 2
     
    dx = (xmax - xmin) / (nx - 1)
    dy = (ymax - ymin) / (ny - 1)
     
    # 初期状態
    p  = np.zeros((ny, nx))
    pd = np.zeros((ny, nx))
    b  = np.zeros((ny, nx))
    c  = np.zeros((ny, nx))
    f  = np.zeros((ny, nx))
    g  = np.zeros((ny, nx))
    h  = np.zeros((ny, nx))
    k  = np.zeros((ny, nx))
    l  = np.zeros((ny, nx))
    o  = np.zeros((ny, nx))
    x  = np.linspace(xmin, xmax, nx)
    y  = np.linspace(xmin, xmax, ny)
     
    # 電荷
    b[ 50  , 1 : 101 : 4 ]  =  -10.e-8
    c[ 50  , 2 : 101 : 8 ]  =  -10.e-8
    f[ 50  , 3 : 101 : 8 ]  =  -10.e-8
    g[ 50  , 4 : 101 : 8 ]  =  -10.e-8
    h[ 47:50 , 1 : 101 : 4 ]  =  -10.e-8
    k[ 47 , 0 : 101 : 8 ]  =  -10.e-8
    l[ 47 , 6 : 101 : 8 ]  =  -10.e-8
    o[ 47 , 7 : 101 : 8 ]  =  -10.e-8
    for it in range(nt):
     
        pd = p.copy()
     
        p[1:-1,1:-1] = (((pd[1:-1, 2:] + pd[1:-1, :-2]) * dy**2 +
                        (pd[2:, 1:-1] + pd[:-2, 1:-1]) * dx**2 -
                        (b[1:-1, 1:-1] + c[1:-1 , 1:-1] + f[1:-1 , 1:-1] + g[1:-1 , 1:-1] + h[1:-1 , 1:-1] + k[1:-1 , 1:-1] + l[1:-1 , 1:-1] + o[1:-1 , 1:-1]) * dx**2 * dy**2) / 
                        (2 * (dx**2 + dy**2)))
     
        p[0, :] = 0
        p[ny-1, :] = 0
        p[:, 0] = 0
        p[:, nx-1] = 0
     
     
    def plot2D(x, y, p):
        fig = plt.figure(figsize=(11, 7), dpi=100)
        ax = fig.gca(projection='3d')
        X, Y = np.meshgrid (x , y)
        surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,linewidth=0, antialiased=False)
        ax.view_init(30, 225)
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        ax.set_zlabel('φ')  
     
    plot2D(x, y, p)
    plt.show()
    
    csv_path = "/content/drive/MyDrive/"
    # CSVファイル名
    csv_name = 'out_name.csv'
    # リストをCSVファイルに出力（'w'の場合は、新規作成）
    with open(csv_path+csv_name, 'w', newline='') as f:
        writer = csv.writer(f)
        writer.writerows([[x, y, p]])   # リスト形式の配列のままcsvファイルに出力
    
    #オブジェクトを生成
    fig = plt.figure(figsize=(11,7), dpi=100)
    fig
     
    xmin = 0
    xmax = 2
    ymin = 0
    ymax = 2
     
    x  = np.linspace(xmin, xmax, nx)
    y  = np.linspace(xmin, xmax, ny)
     
    X, Y = np.meshgrid(x, y)
    plt.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis) 
    plt.colorbar()
    plt.xlabel('X')
    plt.ylabel('Y')
    
[/code]


CC BY-SA 4.0
2022-10-05T11:05:34.383
91458
2022-10-14T03:43:48.547
54717
question
python numpy google-colaboratory
CSVファイルへの出力がうまくいかない


ここにいる人たちはテレパシー初心者なので、なるべくエスパーでなくてもわかるよう、質問者さんのやりたいことを具体的に説明していただければ幸いなのですが、どうも無理そうなので勝手にやりたいことを決めつけた上で回答いたします。ご容赦ください。

2次元格子上での x、y、p(x, y) をそれぞれフラット化して、3列のデータとして出力することにします。

まず、質問者さんが行っているように、2次元[ポアソン方程式](https://ja.wikipedia.org/wiki/%E3%83%9D%E3%82%A2%E3%82%BD%E3%83%B3%E6%96%B9%E7%A8%8B%E5%BC%8F)を[境界値問題](https://ja.wikipedia.org/wiki/%E5%A2%83%E7%95%8C%E5%80%A4%E5%95%8F%E9%A1%8C)として、2次元格子上で差分化し[緩和法](https://en.wikipedia.org/wiki/Relaxation_\(iterative_method\))を用いて解きます。求めるスカラー場が
p(x, y) です（pは、おそらくφの意味）。ポアソン方程式の右辺をb(x, y)とします。境界条件は p(x, y) = 0
です。次のサイトが参考になります。

  * [宇宙に入ったカマキリ - 【第15回Python流体の数値計算】2次元ポアソン方程式をPythonで実装する](https://takun-physics.net/10186/)

[code]

    from matplotlib import pyplot as plt
    import numpy as np
    
    max_iter = 100000  # 最大イテレーション回数
    eps = 1e-4  # イテレーション間の「相対誤差」がこの値以下になったら終了
    
    nx = 101  # x方向の格子点の数
    ny = 101  # y方向の格子点の数
    xmin = 0  # x座標の最小値
    xmax = 2  # x座標の最大値
    ymin = 0  # y座標の最小値
    ymax = 2  # y座標の最大値
    
    dx = (xmax - xmin) / (nx - 1)  # x方向の格子幅
    dy = (ymax - ymin) / (ny - 1)  # y方向の格子幅
    
    p = np.zeros((ny, nx))  # 求めたい場
    b = np.zeros((ny, nx))  # 電荷
    x = np.linspace(xmin, xmax, nx)  # x方向の格子位置
    y = np.linspace(xmin, xmax, ny)  # y方向の格子位置
    
    # 電荷分布を設定する
    
    q0 = -1.0e-7
    
    b[50, 1:101:4] += q0
    b[50, 2:101:8] += q0
    b[50, 3:101:8] += q0
    b[50, 4:101:8] += q0
    b[47:50, 1:101:4] += q0
    b[47, 0:101:8] += q0
    b[47, 6:101:8] += q0
    b[47, 7:101:8] += q0
    
    # 電荷分布をプロットする
    
    im = plt.imshow(b, interpolation="none", extent=(xmin, xmax, ymin, ymax))
    plt.colorbar(im)
    plt.title("charge distribution")
    plt.xlabel("$x$")
    plt.ylabel("$y$")
    plt.show()
    
    # ポアソン方程式を緩和法で解く
    
    for i in range(max_iter):
        p1 = np.zeros((ny, nx))
        p1[1:-1, 1:-1] = (
            (p[1:-1, 2:] + p[1:-1, :-2]) * dy**2
            + (p[2:, 1:-1] + p[:-2, 1:-1]) * dx**2
            - (b[1:-1, 1:-1]) * dx**2 * dy**2
        ) / (2 * (dx**2 + dy**2))
    
        rel_err = np.linalg.norm(p - p1) / np.linalg.norm(p1)
        p = p1
        # if i % 100 == 0:
        #     print(f"i = {i}, rel_err = {rel_err}")
        if rel_err < eps:
            break
    
    print(f"rel_err = {rel_err}")
    
    # 解をプロットする
    
    im = plt.imshow(p, interpolation="none", extent=(xmin, xmax, ymin, ymax))
    plt.colorbar(im)
    plt.title("solution")
    plt.xlabel("$x$")
    plt.ylabel("$y$")
    plt.show()
    
[/code]

![](https://i.stack.imgur.com/ndApb.png)

さて、p(x, y)
は2次元配列として得られていますが、これを[`flatten`メソッド](https://numpy.org/doc/stable/reference/generated/numpy.ndarray.flatten.html)により1次元配列にフラット化して、(x,
y, p(x, y))
の組を1行とする、3列のデータからなる[pandas](https://pandas.pydata.org/)の[データフレーム](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.html)を作成します。

[code]

    import pandas as pd
    
    xx, yy = np.meshgrid(x, y)
    
    df = pd.DataFrame(
        {
            "x": xx.flatten(),
            "y": yy.flatten(),
            "phi": p.flatten(),
        }
    )
    
    print(df)
    
[/code]

[code]

              x    y  phi
    0      0.00  0.0  0.0
    1      0.02  0.0  0.0
    2      0.04  0.0  0.0
    3      0.06  0.0  0.0
    4      0.08  0.0  0.0
    ...     ...  ...  ...
    10196  1.92  2.0  0.0
    10197  1.94  2.0  0.0
    10198  1.96  2.0  0.0
    10199  1.98  2.0  0.0
    10200  2.00  2.0  0.0
    
    [10201 rows x 3 columns]
    
[/code]

あとはデータフレームの[`to_csv`メソッド](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.to_csv.html)を使えばCSVファイルとして保存できます。

[code]

    df.to_csv("output.csv", index=False, header=True)
    
[/code]

お好みで、`index=True`や`header=False`にしてください。

同様に、[`to_excel`メソッド](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.to_excel.html)でExcelファイルとして保存できます（ExcelはCSVを読み込めるので、あまり意味があるかどうかは分かりませんが）。

[code]

    df.to_excel("output.xlsx", index=False, header=True)
    
[/code]


CC BY-SA 4.0
2022-10-14T03:43:48.547
91587
2022-10-14T03:43:48.547
36010
91458
answer