iembry - 1 year ago 242

R Question

I am attempting to rewrite code from GNU Octave/MATLAB in R version 3.3.1. In the original code, A and B were set as global variables in a function and then in the script file both A and B were set as global variables.

In R, this is the error message that I receive when I attempt to use the

`ode45`

Error in eval(expr, envir, enclos) : object 'z1' not found

Can anyone suggest how to set the global variables in R as was done in the GNU Octave/MATLAB code?

Thank you.

R code follows

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

I am using GNU Octave, version 3.8.1 to run the code. The following code in GNU Octave/MATLAB is what I have attempted to replicate above:

`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');

This is

`X`

`X =`

1.0000

1.0016

1.0043

The size of

`t`

`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

The size of

`x`

`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

This is the expected plot:

Recommended for you: Get network issues from **WhatsUp Gold**. **Not end users.**

Answer Source

The following answer uses some elements from the various comments and the answer by 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")
```

Recommended from our users: **Dynamic Network Monitoring from WhatsUp Gold from IPSwitch**. ** Free Download**