# 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. 还有想让两条曲线一条是实线,一条是虚线,也没成功。看论坛里哪位大牛把这两个问题搞一下。
下标的梗,不是bug
调了下画图效果,改成了分析解是线,数值解是点。两个都是线就重合到一起了。
另一个下标的问题好像影响不大?除了画图外,没看见用 t
的值,就没改
# 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 超多示例
好的,我看看哈。非常感谢您!
为什么我执行这个文件时,画不出图,也没任何提示。
但是我在 REPL 下,
julia> x = 0:0.1:1; y = x;
julia> using GR
julia> plot(x,y)
可以正常的画图。
当作脚本运行需要在代码的最后面加上 gui()
或者把 gr()
改成 gr(show = true)
。
我这边是 VSC 可以逐行执行,然后 plot plane 出图,atom 就不行,作为脚本执行也不行,上面那两条也是弹窗,一闪而过。还没有找到好的解决办法,或许可以换个后端。
ref:
好的,谢谢。
现在确实 一闪而过 。。。
实在不行可以存文件嘛,在最后加一句 savefig("my_plot.pdf")
1 个赞
恩,已经可以啦,谢谢