另类方式秀恩爱steemCreated with Sketch.

in #cn7 years ago (edited)

秀恩爱,方式不一而足,不一定要贴照片。比如,有人选择去面馆吃面,有人则要直播吃辣椒。我也凑个热闹。我秀恩爱的方式是——

背景

媳妇甩给我一堆 Matlab 代码:去,给我转成 R 语言。

这是个关于大气化学的常微分方程。Matlab 代码已经解出来了。但是,媳妇考虑到日后维护和拓展的需要,就想换成 R。

娶个女博士做老婆是怎样一种体验?挪,就是这种体验。

好在难不倒我。作为资深 R 语言爱好者,三下五除二,我将它转成了 R 代码。

解是解出来了,但是,得到的解差别很大。

Matlab 提供了一系列解常微分方程的函数,例如 ode15s, ode23s, ode23t, ode23tb 等。我手头代码里用的是 ode23tb。

R 语言的 deSolve 包提供了解常微分方程的函数 ode(),其中可用的参数 method 可以是以下取值:

method = c("lsoda", "lsode", "lsodes", "lsodar", "vode", "daspk",
           "euler", "rk4", "ode23", "ode45", "radau", 
           "bdf", "bdf_d", "adams", "impAdams", "impAdams_d", "iteration")

里面不包括 ode23tb 方法。我写的 R 代码里用了默认的 lsoda 方法,导致得到的解跟 Matlab不同。

我做过的测试

我找了个简单的常微分方程,来自 Matlab 在线帮助

dy/dt = - 10t

我分别用 Matlab 和 R 对它求解。

Matlab 解法:

tspan = [0 2];
y0 = 1;
[t,y] = ode23tb(@(t,y) -10*t, tspan, y0);

R 语言解法:

tspan = seq(0, 2, 0.2)
y0 = 1
testf <- function(t, y, p){
  with(as.list(c(y, p)), {
    dydt = -a*t
    list(dydt)
  })
}
p <- list(a = 10)
zz <- deSolve::ode(y = y0, times = tspan, func = testf, parms = p)

二者得到的结果是不同的。ode() 未指定 method,默认用的是 lsoda。

我尝试了把二者都改为 ode23 方法。即把 Matlab 代码最后一行改为:

[t,y] = ode23(@(t,y) -10*t, tspan, y0);

并把 R 代码最后一行改为:

zz <- deSolve::ode(y = y0, times = tspan, func = testf, parms = p, method = 'ode23')

二者结算的结果就完全一致了。

但是,我手头要改的大气化学常微分方程,可能是因为太复杂,如果用 ode23 方法的话,无论 matlab 还是 R,计算时间都太长,花了12小时都没算完。

我的提问

Matlab 的 ode23tb 函数,在 R 语言里该如何实现?或者反过来,R 语言 deSolve 包 ode() 里的 "lsoda" 方法,如何在 Matlab 里实现?

我知道这个问题我不该问,但是我实在黔驴技穷了,发到这里碰碰运气,看有没有高人指点一下。这个事儿已经连续折磨我两天了,头发掉了一大把。

你说说,这算不算秀恩爱?

Sort:  

有一种网站叫 stackoverflow.

好主意!去那里秀恩爱更另类

你俩属于华山论剑。

重阳一生,不弱于人!

真开心再次见到你!cn区点赞机器人 @cnbuddy 感谢你对cn区作出成长的贡献。 @cnbuddy 旨在助力cn区快速发展,更多cn区动态,请查看我的主页。欢迎关注我们的大股东 @skenan,并注册使用由其开发的 CNsteem.com。倘若你想让我隐形,请回复“取消”。

好奇心殺死貓~

豈止是貓,還有好多腦細胞

数学高考60分的表示祝你好运…^ _ ^

谢谢你的祝福

你可以搞一个类似乌托邦的东西,提问让别人回答,答对了才赞

已经有了类似的东西吧,答对了打赏人民币。改成打赏 SBD 就行了。

可以去乌托邦里提问试试啊。

乌托邦还能提这问题?

有一个Task Requests的选项。。。。 你可以给你的问题包装包装。。。

照你说的,我去问了,看看啥结果。

看看程序员多不多吧,我也有解决不了的问题,准备发上去问问了。。。

这……也真是厉害了……甩代码

比甩硫酸要安全多了

我最喜欢回复这种贴子。
想都不用想。

读都不用读

对啊,但我还是浪费了时间把汉字都读了,好后悔!

可以再浪费一遍,因为我把标题和正文都改了。

你能不能秀一个俺们能看懂的恩爱?

等下回

Coin Marketplace

STEEM 0.17
TRX 0.15
JST 0.028
BTC 62205.55
ETH 2397.85
USDT 1.00
SBD 2.50