[Java] 有限差分法で梁のたわみを求める

# はじめに
Javaと材料力学の復習をかねて作成した。

# 原理
## たわみの基礎式の離散化
分布荷重$p(x)$が加わるヤング率$E$,断面二次モーメント$I$の梁のたわみ$y$は以下の式から求まる。
```math
\frac{d^4y}{dx^4} = \frac{p(x)}{EI}
```

ここで4階の数値微分を以下の式で計算する。
```math
\frac{d^4y}{dx^4} \sim \frac{y_{i-2}-4y_{i-1}+6y_{i}-4y_{i+1}+y_{i+2}}{{\Delta x}^4}
```
よってたわみの基礎式の離散化は以下の式となる。

```math
\frac{y_{i-2}-4y_{i-1}+6y_{i}-4y_{i+1}+y_{i+2}}{{\Delta x}^4} = \frac{p(x_i)}{EI}
```

```math
\frac{1}{{\Delta x}^4}
\begin{bmatrix}
1 & -4 & 6 & -4 & 1\\
\end{bmatrix}
\begin{bmatrix}
y_{i-2} \\
y_{i-1} \\
y_{i} \\
y_{i+1} \\
y_{i+2} \\
\end{bmatrix}
 = \frac{p(x_i)}{EI}
```

節点が$N$個のとき、$i=2$から$i=N-2$までを行列にまとめると以下の式になる。
（※端点は微分の計算で$y_{-2}$や$y_{N+1}$が出てきてしまうので計算できない。）
```math

\frac{1}{{\Delta x}^4}
\begin{bmatrix}
1 & -4 & 6 & -4 & 1 & 0 & &\cdots && 0\\
0 & 1 & -4 & 6 & -4 & 1 & 0 & \cdots && 0\\
 &&&&\vdots &&&&&\\
0 & \cdots && 0 &1 & -4 & 6 & -4 & 1 & 0 \\
0 && \cdots && 0 &1 & -4 & 6 & -4 & 1  \\
\end{bmatrix}
\begin{bmatrix}
y_{1} \\
y_{2} \\
\vdots \\
y_{N-1} \\
y_{N} \\
\end{bmatrix}
 = 
 \frac{1}{EI}
\begin{bmatrix}
p(x_{1}) \\
p(x_{2}) \\
\vdots \\
p(x_{N-1}) \\
p(x_{N}) \\
\end{bmatrix}
```

## 境界条件の離散化
前提として自由端、単純支持、固定端の境界条件は以下の通りとなる。
| 条件 | たわみ$$y$$ | たわみ角$$y'$$ | 曲げモーメント$$EIy''$$ |　せん断応力$$EIy'''$$ |
|:-:|:-:|:-:|:-:|:-:|
| 自由端  |  * |  * |  0 | 0|
| 単純支持  | 0  | *  | 0  |* |
| 固定端  | 0  |  0 | *  | *|

節点$i=0$では前進差分、節点$i=N$では後退差分で離散化する。
以下に節点$i=0$のときの離散化を示す。
節点$i=N$では左辺の行列の非0成分を右端の列にスライドさせればよい。
```math
y(0)=0 \Rightarrow  \begin{bmatrix}
1 &  0 & \cdots & 0\\
\end{bmatrix}
\begin{bmatrix}
y_{1} \\
y_{2} \\
\vdots \\
y_{N-1} \\
y_{N} \\
\end{bmatrix}
 = 0
```
```math
y(0)'=0 \Rightarrow \frac{y_{2}-y_1}{\Delta x} = 0\Rightarrow \begin{bmatrix}
1 & -1 & 0 & \cdots & 0\\
\end{bmatrix}
\begin{bmatrix}
y_{1} \\
y_{2} \\
\vdots \\
y_{N-1} \\
y_{N} \\
\end{bmatrix}
 = 0
```

```math
EIy(0)''=0 \Rightarrow \frac{y_{3}-2y_{2}+y_1}{\Delta x ^2} = 0\Rightarrow \begin{bmatrix}
1 & -2 & 1 & 0 & \cdots & 0\\
\end{bmatrix}
\begin{bmatrix}
y_{1} \\
y_{2} \\
\vdots \\
y_{N-1} \\
y_{N} \\
\end{bmatrix}
 = 0
```

```math
EIy(0)'''=0 \Rightarrow \frac{y_{4}-3y_{3}+3y_{2}-y_1}{\Delta x ^3} = 0\Rightarrow \begin{bmatrix}
1 & -3 & 3 & -1 & 0 & \cdots & 0\\
\end{bmatrix}
\begin{bmatrix}
y_{1} \\
y_{2} \\
\vdots \\
y_{N-1} \\
y_{N} \\
\end{bmatrix}
 = 0
```

## たわみの基礎式+境界条件の行列計算
ここまでの式をまとめて以下の式で表す。
```math
\boldsymbol{A}
\begin{bmatrix}
y_{1} \\
y_{2} \\
\vdots \\
y_{N-1} \\
y_{N} \\
\end{bmatrix}
 = \boldsymbol{b}
```
この式をLU分解を用いた手法で計算したわみ$y_1,y_2,\cdots,y_{N}$ を求める。

# コード
``` Java
package MechDesign.StructureStudy.Beam;

import MechDesign.StructureStudy.StructureStudy;

import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
import org.apache.commons.math4.legacy.linear.ArrayRealVector;
import org.apache.commons.math4.legacy.linear.DecompositionSolver;
import org.apache.commons.math4.legacy.linear.LUDecomposition;
import org.apache.commons.math4.legacy.linear.RealMatrix;
import org.apache.commons.math4.legacy.linear.RealVector;

public class Beam {
    double length;
    int nodeNum;
    double I;
    double E;
    char startBoundaryCondition;
    char endBoundaryCondition;

    public static final class Constraint {
        public static final char FIX = 'F';
        public static final char PIN = 'P';
        public static final char FREE = 'e';
    }

    public Test(double I, double E, double length, int nodeNum, char startBoundaryCondition,
            char endBoundaryCondition) {
        this.I = I;
        this.E = E;
        this.length = length;
        this.nodeNum = nodeNum;
        this.startBoundaryCondition = startBoundaryCondition;
        this.endBoundaryCondition = endBoundaryCondition;
    }

    public double[] solve(double q) {
        // パラメータの設定
        double dx = length / (nodeNum - 1); // メッシュ間隔

        // 係数行列と右辺ベクトルの設定
        RealMatrix A = new Array2DRowRealMatrix(nodeNum, nodeNum);
        RealVector b = new ArrayRealVector(nodeNum);

        // 内部点の設定
        for (int i = 2; i < nodeNum - 2; i++) {
            A.setEntry(i, i - 2, 1 / Math.pow(dx, 4));
            A.setEntry(i, i - 1, -4 / Math.pow(dx, 4));
            A.setEntry(i, i, 6 / Math.pow(dx, 4));
            A.setEntry(i, i + 1, -4 / Math.pow(dx, 4));
            A.setEntry(i, i + 2, 1 / Math.pow(dx, 4));
            b.setEntry(i, q / (E * I));
        }

        // 端点の境界条件
        switch (startBoundaryCondition) {
            case Constraint.FIX:
                A.setSubMatrix(new double[][] { { 1 } }, 0, 0);
                b.setEntry(0, 0.0);
                A.setSubMatrix(new double[][] { { 1 } }, 1, 1);
                b.setEntry(1, 0.0);
                break;
            case Constraint.PIN:
                A.setSubMatrix(new double[][] { { 1 } }, 0, 0);
                b.setEntry(0, 0.0);
                A.setSubMatrix(new double[][] { { 1, -2, 1 } }, 1, 0);
                b.setEntry(1, 0.0);
                break;
            case Constraint.FREE:
                A.setSubMatrix(new double[][] { { 1, -2, 1 } }, 0, 0);
                b.setEntry(0, 0.0);
                A.setSubMatrix(new double[][] { { 1, -3, 3, -1 } }, 1, 0);
                b.setEntry(1, 0.0);
                break;
        }
        switch (endBoundaryCondition) {
            case Constraint.FIX:
                A.setSubMatrix(new double[][] { { 1 } }, nodeNum - 2, nodeNum - 1);
                b.setEntry(nodeNum - 2, 0.0);
                A.setSubMatrix(new double[][] { { 1 } }, nodeNum - 1, nodeNum - 2);
                b.setEntry(nodeNum - 1, 0.0);
                break;
            case Constraint.PIN:
                A.setSubMatrix(new double[][] { { 1 } }, nodeNum - 2, nodeNum - 1);
                b.setEntry(nodeNum - 2, 0.0);
                A.setSubMatrix(new double[][] { { 1, -2, 1 } }, nodeNum - 1, nodeNum - 3);
                b.setEntry(nodeNum - 1, 0.0);
                break;
            case Constraint.FREE:
                A.setSubMatrix(new double[][] { { 1, -2, 1 } }, nodeNum - 2, nodeNum - 3);
                b.setEntry(nodeNum - 2, 0.0);
                A.setSubMatrix(new double[][] { { 1, -3, 3, -1 } }, nodeNum - 1, nodeNum - 4);
                b.setEntry(nodeNum - 1, 0.0);
                break;
        }
        // 連立方程式を解く
        DecompositionSolver solver = new LUDecomposition(A).getSolver();
        RealVector solution = solver.solve(b);

        return solution.toArray();
    }

```

```Java
    public static void main(String[] args) {
        Beam beam = new Beam(1e-4, 10., 50, Beam.Constraint.FIX, Beam.Constraint.FIX, 210e-9);
        double[] yArr = beam.solve(1000);
        double max = yArr[0];
        for (var y : yArr) {
            System.out.println(y * 1e3);
            if (y > max)
                max = y;
        }

        System.out.println("max:" + max);
    }
```