用julia写的用有限元求解1DBVP

综合讨论区
1 
# 1D有限元求解示例。(在0.6与1.0下均可运行)
# 节点编号为：0, 1, ..., n+1
# 0与n+1为边值，待求量只有1, 2, ..., n
# 所求BVP
# y'' = 4y
# y(0) = 1.
# y(1) = e^2

using Plots
gr()

function BVPFEM(inter, bv, n)  # inter为计算区间; bv为边值; n为步数

    a1 = inter[1]; a2 = inter[2]
    ya = bv[1];    yb = bv[2]
    h  = (a2-a1) / (n+1)

    alpha = (8. / 3.)*h + 2. / h
    beta  = 2. * h / 3. - 1. / h

    # 左端项系数矩
    A = zeros(n,n)
    for i in 1:n
        A[i,i] = alpha
    end
    for i in 2:n
        A[i-1,i] = beta
        A[i,i-1] = beta
    end

    # 变量t
    t = zeros(n+2)
    t = [a1+(i-1.)*h for i in 1:n+2]

    # 右端项
    B = zeros(n,1)
    B[1,1] = -ya*beta
    B[n,1] = -yb*beta
    c = zeros(n+2,1)
    c[1,1] = ya; c[n+2,1] = yb
    c[2:n+1,1] = A\B
    # println(c)  # 原方程的数值解
    c1 = exp.(2*t)  # 原方程的解析解
    # println(c1)

    return t, c, c1
end

n = 99  # n为待求量个数
t, c, c1 = BVPFEM([0. 1.], [1. exp.(2.0)], n)  # [0. 1.]为计算区间; [1. exp.(2.0)]为边界值

plot(t, c)
plot!(t, c1)

# 程序有点小瑕疵。1. 数组t的下标为1:n+2，本来想变化到0:n+1(为了与理论推导时的下表一致)，却失败了。2. 还有想让两条曲线一条是实线，一条是虚线，也没成功。看论坛里哪位大牛把这两个问题搞一下。
2

下标的梗，不是bug

3

调了下画图效果，改成了分析解是线，数值解是点。两个都是线就重合到一起了。
另一个下标的问题好像影响不大？除了画图外，没看见用 t 的值，就没改

1DBVP
1DBVP1105×742 77 KB

# 1D有限元求解示例。(在0.6与1.0下均可运行)
# 节点编号为：0, 1, ..., n+1
# 0与n+1为边值，待求量只有1, 2, ..., n
# 所求BVP
# y'' = 4y
# y(0) = 1.
# y(1) = e^2

using Plots
gr()

function BVPFEM(inter, bv, n)  # inter为计算区间; bv为边值; n为步数

    a1 = inter[1]; a2 = inter[2]
    ya = bv[1];    yb = bv[2]
    h  = (a2-a1) / (n+1)

    alpha = (8. / 3.)*h + 2. / h
    beta  = 2. * h / 3. - 1. / h

    # 左端项系数矩
    A = zeros(n,n)
    for i in 1:n
        A[i,i] = alpha
    end
    for i in 2:n
        A[i-1,i] = beta
        A[i,i-1] = beta
    end

    # 变量t
    t = zeros(n+2)
    t = [a1+(i-1.)*h for i in 1:n+2]

    # 右端项
    B = zeros(n,1)
    B[1,1] = -ya*beta
    B[n,1] = -yb*beta
    c = zeros(n+2,1)
    c[1,1] = ya; c[n+2,1] = yb
    c[2:n+1,1] = A\B
    # println(c)  # 原方程的数值解
    c1 = exp.(2*t)  # 原方程的解析解
    # println(c1)

    return t, c, c1
end

n = 99  # n为待求量个数
t, c, c1 = BVPFEM([0. 1.], [1. exp.(2.0)], n)  # [0. 1.]为计算区间; [1. exp.(2.0)]为边界值

plot(t, c,
        markershape = :circle,
        markersize = 3,
        label="Numerical sol", legend = :topleft)
plot!(t, c1,
        linewidth = 2,
        label="Analytical sol")

# 程序有点小瑕疵。
# 1. 数组t的下标为1:n+2，本来想变化到0:n+1(为了与理论推导时的下表一致)，却失败了。
# 2. 还有想让两条曲线一条是实线，一条是虚线，也没成功。看论坛里哪位大牛把这两个问题搞一下。

画图常用ref，前两个是 plots 的文档，后一个是个 gist 超多示例

4 
好的，我看看哈。非常感谢您！
5

为什么我执行这个文件时，画不出图，也没任何提示。
但是我在 REPL 下，

julia> x = 0:0.1:1; y = x;
julia> using GR
julia> plot(x,y)

可以正常的画图。

6

当作脚本运行需要在代码的最后面加上 gui() 或者把 gr() 改成 gr(show = true)。

我这边是 VSC 可以逐行执行，然后 plot plane 出图，atom 就不行，作为脚本执行也不行，上面那两条也是弹窗，一闪而过。还没有找到好的解决办法，或许可以换个后端。

ref:

7

好的，谢谢。
现在确实 一闪而过 。。。

8

实在不行可以存文件嘛，在最后加一句 savefig("my_plot.pdf")

1 个赞
9

恩，已经可以啦，谢谢 :grin: