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

4 次查看（过去 30 天）
subbu
subbu 2021-1-23
编辑： Stephen23 2021-1-27
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)];

请先登录，再回答此问题。

采纳的回答

Stephen23
Stephen23 2021-1-27
编辑：Stephen23 2021-1-27
Ran in:
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:
https://www.mathworks.com/help/matlab/math/parameterizing-functions.html
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

更多回答（1 个）

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

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by