bsm公式的matlab代码_Matlab真的不可替代么 迁移Octave实战

说在前面的话

光说不练假把式,这次来点干货,我们实战一次。

7年前还在读博士的时候,在科学网写过一篇博客,用Matlab以不超过150行代码实现了一个简单的二维弹性有限元程序,我们不妨拿这个小例子开刀(Matlab代码可以从 http://blog.sciencenet.cn/blog-3760-678647.html 下载)。

我这个小程序没法在Octave下直接运行,因为涉及到了PDE Toolbox的使用。很多人说Matlab不可替代的是Toolbox工具箱,我也想在这篇文章argue一下,不是所有的工具箱都不能替代的,其实作为Matlab的替代品,Octave也仿照Matlab构建了大量的package,大部分是第三方作者开源贡献的,不信您上 https://octave.sourceforge.io/packages.php 瞧瞧,是不是有很多相似的面孔。我想借用毛主席的名言,“没有调查就没有发言权” (:

Matlab PDE工具箱的替代

我们今天为了完成这个小迁移就来找找PDE Toolbox的替代品。

Octave用来替代PDE Toolbox的工具包主要有两个,一个是生成网格用的msh(https://octave.sourceforge.io/msh/index.html),一个是显示网格和计算结果用的fpl(https://octave.sourceforge.io/fpl/index.html)。

先来看看生成网格的msh。

msh生成非结构化网格实际是后台调用的gmsh,刚好我的macbook上已经有了,如果没有该程序可以简单的用homebrew安装一下,执行:

brew install gmsh

安装工具包的过程很简单,将octave目录转移到下载的工具包所在目录后,执行pkg install xxx,我们这儿执行:

pkg install splines-1.3.3.tar.gz

另外msh工具包还要依赖于spline样条函数工具包,也需要提前下载好,同样执行:

pkg install msh-1.0.10.tar

fpl工具包安装也类似,同样执行:

pkg install fpl-1.3.5.tar.gz

至此,完成我们所有准备工作。

Matlab程序迁移

一个不幸的消息是Matlab工具包PDE Toolbox中的initmesh函数没有直接的替代品,需要使用msh工具包中的函数改写。

fpl的官方文档即给出了和我们Matlab程序生成三角形网格类似的例子,直接抄过来,如下图所示:

6c7df100ec25b8b49fbc818d66ca73d8.png

Octave改写的网格生成程序

熟悉gmsh的同学可以发现,这个程序就是简单的调用了一下gmsh而已,连命令都一样(: 

执行该程序可以得到我们之前在Matlab中一样的网格,数据结构也是和Matlab PDE工具箱中的一样,结果如下图所示:

a21c846f7721f134708f1d64a112f403.png

Octave生成的二维网格

好了,Matlab中的initmesh函数用我们上面写的程序换换就好。plf工具包生成的网格都在mesh结构体变量里,我们把它简单的赋值出来给p、t、e就好,这也是我们这个小程序移植过程中唯一的改动。

除此以外,还有个需要注意的地方,Octave里第三方工具箱使用前需要pkg load一下。

好了,我们把所有的改动都写在一起:

6e150922bd83a9f115b3666c3cc631be.png

Matlab程序的改写

运行一下试试,然后我们画图:

8f2e74acdc7ce8a0e955e42effa252fb.png Octave画图代码

159831c9b578ae089f8b65a5abba922f.png

Octave画图结果

这就完成了Matlab程序的迁移了,是不是很简单。

所以说,Matlab被禁大家不要怕,迁移很容易,如果您不会,那就找我呀(:

用FELAC写一个二维弹性有限元程序

最后我还想说,其实写一个弹性力学有限元程序还有更简单的方法 - 

这就是用我们有限元语言及其编译器FELAC软件。

管你MATLAB是不是被禁,管你OCTAVE是不是有相关的工具箱,用FELAC写有限元程序真的只要不到一分钟 ^-^

bef9d0288512ebd718214de3760daac4.png

FELAC公式库里的FELAC二维线弹性有限元语言脚本

d7058830b43519f97a2fec1e0050865c.png

FELAC生成的二维线弹性有限元C语言主程序

db4925ee36ff3402a6243f936878788f.png

FELAC计算结果,x方向位移

9ffe624500881eb4ebd026a4574cc7db.png

FELAC计算结果,y方向位移

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值