R - 使用全局變量來解決類似於MATLAB/GNU的ode45 Octave

Question

我試圖在R版本3.3.1中重寫來自GNU Octave/MATLAB的代碼。在原始代碼中，A和B在函數中被設置爲全局變量，然後在腳本文件中A和B都被設置爲全局變量。

在R，這是我接收時，我試圖使用ode45函數的錯誤消息：

錯誤的eval（expr中，ENVIR，enclos）：對象 'Z1' 沒有找到

任何人都可以提出如何在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八度，版本3.8.1運行代碼。在GNU八度/ 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 

這是預期的情節：

Answer 1

以下答案使用來自各種評論和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") 

Answer 2

我相信這是你想要的東西：

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") 

所做的更改：

1. 使用矩陣乘法運算符在GRT，而不是標量乘法
2. 固定的X0索引（稱爲x在你的函數）
3. 註釋掉2條不必要從GRT
4. 增加了return陳述線GRT

我不得不做出一些假設您在語法錯誤的地方嘗試做什麼，如索引X0。由於我沒有octave引用的示例輸出圖（並且我無法讓您的代碼在我的Octave CLI中運行），所以我無法判斷這些假設是否正確，如果不是，我的情節可能會有所不同。

這是從上面的代碼生成的情節：

最後一點：看起來你從來沒有使用過的P <- mat("t, x")結果無論如何，但我不認爲它做什麼，你認爲這是無論如何，基於結果的對象。