R - 使用全局变量求解 ode45 类似于 MATLAB/GNU Octave
R - use global variables to solve ode45 similar to MATLAB/GNU Octave
我正在尝试在 R 版本 3.3.1 中重写来自 GNU Octave/MATLAB 的代码。原代码中在函数中将A和B设置为全局变量，然后在脚本文件中将A和B都设置为全局变量
在 R 中，这是我尝试使用 ode45 函数时收到的错误消息：
Error in eval(expr, envir, enclos) : object 'z1' not found
谁能建议如何像在 GNU Octave/MATLAB 代码中那样在 R 中设置全局变量？
谢谢。
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：
X =
1.0000
1.0016
1.0043
t 的大小是 809 行，1 列。这是 t 的部分内容。
t =
0.00000
0.00305
0.00632
0.00928
0.01226
0.01524
0.01840
0.02186
0.02482
0.02778
0.03079
0.03391
0.03750
0.04046
0.04344
0.04646
0.04959
0.05321
0.05618
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
这是预期的情节：
我相信这就是你想要的：
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 条不必要的行
- 向
GRT 添加了 return 语句
我不得不对你在语法错误的地方尝试做的事情做出一些假设，比如 X0 的索引。由于我没有来自 Octave 的示例输出图可供参考（而且我无法在我的 Octave CLI 中将您的代码获取到 运行）我无法判断这些假设是否正确，如果不是我的情节可能不同。
这是上面代码的结果图：
最后的说明：看起来你从来没有使用过 P <- mat("t, x") 的结果，但我不认为它会根据结果对象做你认为它正在做的事情。
以下答案使用了各种评论中的一些元素和 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")
我正在尝试在 R 版本 3.3.1 中重写来自 GNU Octave/MATLAB 的代码。原代码中在函数中将A和B设置为全局变量，然后在脚本文件中将A和B都设置为全局变量
在 R 中，这是我尝试使用 ode45 函数时收到的错误消息：
Error in eval(expr, envir, enclos) : object 'z1' not found
谁能建议如何像在 GNU Octave/MATLAB 代码中那样在 R 中设置全局变量？
谢谢。
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：
X =
1.0000
1.0016
1.0043
t 的大小是 809 行，1 列。这是 t 的部分内容。
t =
0.00000
0.00305
0.00632
0.00928
0.01226
0.01524
0.01840
0.02186
0.02482
0.02778
0.03079
0.03391
0.03750
0.04046
0.04344
0.04646
0.04959
0.05321
0.05618
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
这是预期的情节：
我相信这就是你想要的：
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 条不必要的行
- 向
GRT 添加了
return 语句
我不得不对你在语法错误的地方尝试做的事情做出一些假设，比如 X0 的索引。由于我没有来自 Octave 的示例输出图可供参考（而且我无法在我的 Octave CLI 中将您的代码获取到 运行）我无法判断这些假设是否正确，如果不是我的情节可能不同。
这是上面代码的结果图：
最后的说明：看起来你从来没有使用过 P <- mat("t, x") 的结果，但我不认为它会根据结果对象做你认为它正在做的事情。
以下答案使用了各种评论中的一些元素和 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")