R - 使用全局变量来解决类似于 MATLAB/GNU Octave 的 ode45答案

【问题标题】：R - use global variables to solve ode45 similar to MATLAB/GNU OctaveR - 使用全局变量来解决类似于 MATLAB/GNU Octave 的 ode45
【发布时间】：2016-07-30 15:53:29
【问题描述】：

我正在尝试从 R 版本 3.3.1 中的 GNU Octave/MATLAB 重写代码。在原始代码中，A 和 B 被设置为函数中的全局变量，然后在脚本文件中，A 和 B 都被设置为全局变量。

在 R 中，这是我尝试使用 ode45 函数时收到的错误消息：

eval(expr, envir, enclos) 中的错误：找不到对象“z1”

谁能建议如何在 R 中设置全局变量，就像在 GNU Octave/MATLAB 代码中所做的那样？

谢谢。

R 代码如下

#list.R is a wrapper for list used to replicate the GNU Octave/MATLAB syntax


source("https://raw.githubusercontent.com/ggrothendieck/gsubfn/master/R/list.R")

install.load::load_package("ramify", "pracma")

GRT <- function (t, x) {
         A  <- A
         B  <- B
         z1 <- x[1, 1]; z2 <- x[2, 1]; z3 <- x[3, 1]; X <- mat("z1; z2; z3")
         xd <- A * X + B * exp(-t) * sin(400 * t)
}


A <- -mat("2, 3, 2; 1, 5, 3; 2, 3, 1")
B <- mat("1; 3; 2")

ts <- 0.0
tf <- 10
X0 <- c(1, 0, -1)

list[t, x] <- ode45(GRT, ts, tf, X0)

P <- mat("t, x")

matplot(t, x, xlab = "time - (s)", ylab = "x")

我正在使用 GNU Octave 版本 3.8.1 来运行代码。 GNU Octave/MATLAB 中的以下代码是我在上面尝试复制的：

function xd=GRT(t,x)
global A B
z1=x(1,1); z2=x(2,1);z3=x(3,1);X=[z1;z2;z3];
xd =A*X+B*exp(-t)*sin(400*t);
endfunction % only needed for GNU Octave

global A B
A = -[2,3,2;1,5,3;2,3,1];
B = [1;3;2];
ts = 0.0;
tf = 10;
T=[ts,tf]; X0=[1,0,-1];
[t,x] = ode45(@GRT,T,X0)
P = [t,x];
plot(t,x)
xlabel('time - (s)');
ylabel('x');

这是X：

t 的大小为 809 行 1 列。这是对t 的部分介绍。

x 的大小为 809 行，3 列。这是对x 的部分介绍。

x =

1.0000e+00   0.0000e+00  -1.0000e+00
1.0016e+00   1.0937e-02  -9.9982e-01
1.0043e+00   2.5810e-02  -9.9752e-01
1.0040e+00   3.1460e-02  -1.0007e+00
1.0012e+00   2.9337e-02  -1.0090e+00
9.9908e-01   2.9132e-02  -1.0161e+00
1.0001e+00   3.8823e-02  -1.0170e+00
1.0028e+00   5.4307e-02  -1.0148e+00
1.0026e+00   6.0219e-02  -1.0178e+00
9.9979e-01   5.8198e-02  -1.0260e+00
9.9739e-01   5.7425e-02  -1.0336e+00
9.9809e-01   6.6159e-02  -1.0351e+00
1.0007e+00   8.1786e-02  -1.0331e+00
1.0004e+00   8.7669e-02  -1.0361e+00
9.9753e-01   8.5608e-02  -1.0444e+00
9.9500e-01   8.4599e-02  -1.0522e+00

这是预期的情节：

【问题讨论】：

首先，警告——在编程中使用全局变量被认为是非常糟糕的做法。话虽如此，您能否展示您的预期结果（Octave 中的情节）？
我不熟悉 R，但在 MATLAB 中 @Hack-R 所说的非常正确。它非常容易出错并且非常难以调试（例如，您设置了全局值，然后在其他一些函数中您不小心覆盖了它，然后第一个函数出错了，但不知道为什么。）
所以在那个函数GRT 中你使用的是x，我认为它是X0，但是你试图通过索引引用的维度没有意义，因为@987654337 @ 只是一个 3 元素向量。
还有你想将哪个对象作为t传递给GRT？ ts 或 tf?
似乎mat() 在使用函数调用时无法找到变量z1 等。尝试改用X <- matrix(c(z1, z2, z3), byrow = TRUE, nrow = 3)。然后我假设A * X 应该是矩阵乘法。这在 R 中写为A %*% X。如果我进行这两项更改，代码就会运行。输出是否是预期的，我不能说...

标签： r matlab global-variables octave

【解决方案1】：

我相信这就是你想要的：

source("https://raw.githubusercontent.com/ggrothendieck/gsubfn/master/R/list.R")

pacman::p_load(ramify, pracma) # I use pacman, you don't have to

GRT <- function (t, x) {
  X <- mat("z1; z2; z3")
  xd <- A %*% X + B %*% exp(-t) * sin(400 * t)
  return(z1)
}


A <- -mat("2, 3, 2; 1, 5, 3; 2, 3, 1")
B <- mat("1; 3; 2")

ts <- 0.0
tf <- 10
X0 <- c(1, 0, -1)
z1 <- X0[1]
z2 <- X0[2]
z3 <- X0[3] 

GRT(t=ts,x=X0)

list[t, x] <- ode45(GRT, ts, tf, X0)

P <- mat("t, x")

matplot(t, x, xlab = "time - (s)", ylab = "x")

所做的更改：

在GRT 中使用矩阵乘法运算符而不是标量乘法
修复了 X0 的索引（在您的函数中称为 x）
从GRT 中注释掉了 2 行不必要的行
将return 语句添加到GRT

我不得不对你在语法错误的地方尝试做什么做出一些假设，比如X0 的索引。由于我没有 Octave 的示例输出图可供参考（而且我无法让您的代码在我的 Octave CLI 中运行），我无法判断这些假设是否正确，如果不是，我的图可能会有所不同。

这是上面代码的结果图：

最后说明：您似乎从未使用过P <- mat("t, x") 的结果，但我认为它并没有按照您认为的那样做，基于结果对象。

【讨论】：

【解决方案2】：

以下答案使用了来自各种 cmets 的一些元素以及 Hack-R 的答案。

source("https://raw.githubusercontent.com/ggrothendieck/gsubfn/master/R/list.R")

install.load::load_package("ramify", "pracma")

GRT <- function (t, x) {

z1 <- x[1, 1]; z2 <- x[2, 1]; z3 <- x[3, 1]; X <- matrix(data = c(z1,  
z2, z3), nrow = 3, ncol = 1)

xd <- A %*% X + B * exp(-t) * sin(400 * t)

return(xd)
}


A <- -mat("2, 3, 2; 1, 5, 3; 2, 3, 1")

B <- mat("1; 3; 2")

ts <- 0.0

tf <- 10

X0 <- c(1, 0, -1)

list[t, x] <- ode45(GRT, ts, tf, X0, atol = 0.000001, hmax = 1.0)

matplot(t, x, xlab = "time - (s)", ylab = "x", type = "l")

【讨论】：