# how to extract an extra variable which is calculated in ode function to main workspace?

1 view (last 30 days)
subbu on 23 Jan 2021
Edited: Stephen on 27 Jan 2021
I have to extract the values of "p" variable which is calculated in 'model' function. can anyone help?
clear;
clc;
ifp = fopen('loma.txt','r');
% capture and discard the first two lines of the file
for i=1:46
R = fgetl(ifp); % delete redundant lines
end
% load the acceleration data 90 degree
data90 = fscanf(ifp,'%f')/100; % m/s/s
% close the input file
fclose(ifp);
ifp = fopen('loma.txt','r');
% capture and discard the first two lines of the file
for i=1:1646
F = fgetl(ifp); % delete redundant lines
end
% load the acceleration data 0 degree
data0 = fscanf(ifp,'%f')/100; % m/s/s
% close the input file
fclose(ifp);
dt = 0.02;
npoint=length(data90);
time=0:dt:(npoint-1)*dt;
t=time';
kx=11*10^4; %N/m
ky=8*10^4; %N/m
fyxp=550*10^0; %N
fyxn=650*10^0; %N
fyyp=440*10^0; %N
fyyn=520*10^0; %N
alphax=0.01;
alphay=0.01;
zux0p=fyxp/kx;
zuy0p=fyyp/ky;
zux0n=fyxn/kx;
zuy0n=fyyn/ky;
Z0(1)=0;
Z0(2)=0;
Z0(3)=0;
Z0(4)=0;
Z0(5)=0;
Z0(6)=0;
[t,Z]=ode23s(@model,t,Z0,[],data0,data90,t);
qx=alphax*kx*Z(:,1)+(1-alphax)*ky*Z(:,5);
qy=alphay*ky*Z(:,3)+(1-alphay)*zuy0p/zux0p*ky*Z(:,6);
function [ydot]=model(t,Y,ax2,ay2,t1)
ax = interp1(t1,ax2,t);
ay = interp1(t1,ay2,t);
kx=11*10^4; %N/m
ky=8*10^4; %N/m
fyxp=550*10^0; %N
fyxn=650*10^0; %N
fyyp=440*10^0; %N
fyyn=520*10^0; %N
alphax=0.01;
alphay=0.01;
beta=0.7;
dt = 0.02;
eta=1;
nu=1;
m=50;
xi=5/100;
s=8;
sigma=.2;
mu=0.1;
zux0p=fyxp/kx;
zuy0p=fyyp/ky;
zux0n=fyxn/kx;
zuy0n=fyyn/ky;
zu=((1+sign(Y(5)))*abs(zux0p)+(1-sign(Y(5)))*abs(zux0n))/2/nu;
Cx=Y(5)/zu*(1+beta*(sign(Y(2)*Y(5))-1));
Cy=Y(6)/zu*(1+beta*(sign(Y(4)*Y(6))-1));
px=sqrt(2/pi)*s/sigma*exp(-0.5*((sign(Y(2))*Y(5)/zu-mu)/sigma)^2);
py=sqrt(2/pi)*s/sigma*exp(-0.5*((sign(Y(4))*Y(6)/zu-mu)/sigma)^2);
H11=eta+(1-Y(5)/zu*Cx)*px;
H12=-Y(5)/zu*Cy*py;
H21=-Y(6)/zu*Cx*px;
H22=eta+(1-Y(6)/zu*Cy)*py;
G1=(1-Y(5)/zu*Cx)*Y(2)-Y(5)*zux0p/zu/zuy0p*Cy*Y(4);
G2=-Y(6)/zu*Cx*Y(2)+(1-Y(6)/zu*Cy)*zux0p/zuy0p*Y(4);
Ydot(5)=1/(H11*H22-H12*H21)*(H22*G1-H12*G2);
zdashx=1/(H11*H22-H12*H21)*(H22*G1-H12*G2);
Ydot(6)=1/(H11*H22-H12*H21)*(H11*G2-H21*G1);
zdashy=1/(H11*H22-H12*H21)*(H11*G2-H21*G1);
Ydot(1)=Y(2);
Ydot(3)=Y(4);
qx=alphax*kx*Y(1)+(1-alphax)*ky*Y(5);
qy=alphay*ky*Y(3)+(1-alphay)*zuy0p/zux0p*ky*Y(6);
Ydot(2)=-ax-2*xi*sqrt(kx/m)*Y(2)-qx/m;
Ydot(4)=-ay-2*xi*sqrt(ky/m)*Y(2)-qy/m;
p=1/zu^2*(Y(5)*(Y(2)-mu*zdashx)+Y(6).*(zux0p/zuy0p*Y(4)-mu*zdashy))
% v=cumtrapz(dt,p);
% Energy=(1-alphax)/2.*v;
ydot=[Ydot(1);Ydot(2);Ydot(3);Ydot(4);Ydot(5);Ydot(6)];

Stephen on 27 Jan 2021
Edited: Stephen on 27 Jan 2021
Ignore advice that you should use assignin or evalin, because this does not take into account how ODE solvers work, as well as being slow and complex (this topic has been discussed many times on this forum).
Here is a simple way to retrieve the p value from your function for exactly the t and Z values returned by the ODE solver.
Unlike the incomplete advice to use assignin/evalin here is a complete, simple, efficient, working solution using exactly the method that I gave you in my earlier comment. I also changed your obsolete syntax for providing extra inputs to the function with the recommended function parameterization:
I marked the altered lines with !!! so you can find them easily:
ifp = fopen('loma.txt','r');
% capture and discard the first two lines of the file
for i=1:46
fgetl(ifp); % delete redundant lines
end
% load the acceleration data 90 degree
data90 = fscanf(ifp,'%f')/100; % m/s/s
% close the input file
fclose(ifp);
ifp = fopen('loma.txt','r');
% capture and discard the first two lines of the file
for i=1:1646
fgetl(ifp); % delete redundant lines
end
% load the acceleration data 0 degree
data0 = fscanf(ifp,'%f')/100; % m/s/s
% close the input file
fclose(ifp);
dt = 0.02;
npoint=length(data90);
time=0:dt:(npoint-1)*dt;
kx=11*10^4; %N/m
ky=8*10^4; %N/m
fyxp=550*10^0; %N
fyxn=650*10^0; %N
fyyp=440*10^0; %N
fyyn=520*10^0; %N
alphax=0.01;
alphay=0.01;
zux0p=fyxp/kx;
zuy0p=fyyp/ky;
zux0n=fyxn/kx;
zuy0n=fyyn/ky;
Z0 = [0;0;0;0;0;0];
fun = @(T,X)model(T,X,data0,data90,time(:));% !!! replace obsolete syntax !!!
[t_out,Z_out] = ode23s(fun,time(:),Z0)% !!! and use different variable names !!!
t_out = 2000×1
0 0.0200 0.0400 0.0600 0.0800 0.1000 0.1200 0.1400 0.1600 0.1800
Z_out = 2000×6
0 0 0 0 0 0 -0.0000 -0.0010 0.0000 0.0001 -0.0000 0.0000 -0.0000 -0.0010 -0.0000 -0.0005 -0.0000 -0.0000 -0.0000 0.0002 -0.0000 -0.0011 -0.0000 -0.0000 -0.0000 0.0022 -0.0000 -0.0015 -0.0000 -0.0000 0.0000 0.0044 -0.0001 -0.0019 0.0000 -0.0000 0.0002 0.0061 -0.0001 -0.0030 0.0000 -0.0000 0.0003 0.0065 -0.0002 -0.0054 0.0000 -0.0000 0.0004 0.0054 -0.0003 -0.0083 0.0000 -0.0000 0.0005 0.0032 -0.0005 -0.0104 0.0000 -0.0000
[~,p] = cellfun(fun, num2cell(t_out), num2cell(Z_out,2), 'uni',0); % !!! p values
p = cell2mat(p) % !!! and here is your output p vector !!!
p = 2000×1
0 0.0000 0.0000 0.0000 0.0000 0.0004 0.0017 0.0038 0.0061 0.0083
subbu on 27 Jan 2021
Thanks for the help. I really appreciate it

Sai Veeramachaneni on 27 Jan 2021
Hi,
You can leverage assignin and evalin functions for your usecase.
Hope it helps.

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!