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

7 ビュー (過去 30 日間)
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 日
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

その他の回答 (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 ExchangeMathematics についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by