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

Question

subbu 2021-1-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/725087-how-to-extract-an-extra-variable-which-is-calculated-in-ode-function-to-main-workspace

编辑： 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)];

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Stephen23 2021-1-23

Just return it as an output argument, as shown here:

https://www.mathworks.com/matlabcentral/answers/512891-passing-extra-output-arguments-from-ode

https://www.mathworks.com/matlabcentral/answers/429179-how-do-i-extract-an-intermediate-variable-calculated-and-used-inside-my-ode45-function

Avoid obsolete, undocumented syntaxes by parameterizing the function properly:

https://www.mathworks.com/matlabcentral/answers/517469-some-quetions-about-ode45#answer_426013

https://www.mathworks.com/matlabcentral/answers/469314-ode45-varargin-additional-arguments-after-options#answer_381161

The recommended approach (e.g. using an anonymous function):

https://www.mathworks.com/help/matlab/math/parameterizing-functions.html

请先登录，再进行评论。

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

Answer 1

Stephen23 2021-1-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/725087-how-to-extract-an-extra-variable-which-is-calculated-in-ode-function-to-main-workspace#answer_607268

编辑：Stephen23 2021-1-27

在 MATLAB Online 中打开

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