首页 文章

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

提问于
浏览
0

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

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

eval(expr,envir,enclos)中的错误:找不到对象'z1'

任何人都可以建议如何在GN中设置全局变量,如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

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

这是预期的情节:

GRT plot

2 回答

  • 0

    我相信这就是你想要的:

    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 CLI中运行,我无法判断这些假设是否正确,如果不是,我的情节可能会有所不同 .

    这是上面代码的结果图:

    enter image description here

    最后的注意事项:看起来你从来没有使用过 P <- mat("t, x") 的结果,但我还是根据产生的对象做了什么 .

  • 3

    以下答案使用了各种评论中的一些元素以及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")
    

    ode45 plot in R

相关问题