How to (efficiently) replace eval() with subs() for symbolic equations

Hi,

I've been trying to convert a code I have for modeling a curve into a Simulink model (so, eventually, I'll be able to run automatic parameter estimation). The code itself is comprised of a few different symbolic functions. I previously used the code in Matlab 2012b (where it exectutes very quickly without problems) however I couldn't make it work in Simulink because of the C++ compiler issue. I tried to apply the patches to resolve the compiler problem with no sucess so I gave up and switced over to 2018b. Now the problem seems to be trying to replace the eval() function which is no longer supported. My best solution to date uses sub() and str2sym(). This works but is so slow that it's effectively useless. Does anybody have suggestions on how I can improve the code?

Old version:

for kt=1:nnt       
   tt=radt(kt); s=alfa/tt;  bt=beta/tt; btF=bt.*eval(FF);  ft(kt)=sum(real(btF));   
end

New version:

btf=@(s,bt) double(bt.*subs(str2sym(Fs))); 
for kt=1:nnt                  
    tt=radt(kt); s=alfa/tt; bt=beta/tt; btF=btf(s,bt); ft(kt)=sum(real(btF)); 
end

For reference, a simplified version of FF is: '(0.1*s)^(1/2)*besselk(0.0, 8.0*(0.1*s)^(1/2))'

I'll include the files being used in case anyone would like to take a look in more detail.

Thanks in advance for any assistance!

[Walter Roberson]

 btf = matlabFunction(bt*str2sym(Fs), 'vars', {s, bt})

[Nat A]

Thanks for the help, much appreciated.

When I replace the first line in the new version of the code with the one suggested above I get the following error:

Undefined function or variable 's'.
Error in INVLAP_2018_b (line 31)
 btf = matlabFunction(bt*str2sym(Fs), 'vars', {s, bt})
Error in N_89_b (line 24)
[t,hob]=INVLAP_2018_b(h_ob_b,tmin,tmax,tnum,6,20,19);

The s I'm using is a Laplace variable so it was never specifically defined in the code. I went back and tried adding

syms s

in the function file (h_ob_b.m) that is called when assigning Fs but I ended up with the same error message. The full function assigned to Fs is:

function h_ob_b = h_OB_b(S_est,T_est,Co,Cs,rw,lamda,r,Ho)
phi = sprintf('%e./%e.*s',S_est,T_est);
E_3 = sprintf('(%e.*s.*besselk(0,sqrt(%s).*%e)+%e.*sqrt(%s).*besselk(1,sqrt(%s).*%e))./(%e.*s)',Co,phi,rw,lamda,phi,phi,rw,Co);
E_4 = sprintf('(%e.*s.*besselk(0,sqrt(%s).*%e)+%e.*sqrt(%s).*besselk(1,sqrt(%s).*%e))./(%e.*s)',Cs,phi,rw,lamda,phi,phi,rw,Cs);
gamma_2 = sprintf('besselk(0,sqrt(%s).*%e).^2-(%s).*(%s)',phi,r,E_3,E_4);
h_ob_b = sprintf('-%e.*%e.*sqrt(%s).*besselk(0,sqrt(%s).*%e)*besselk(1,sqrt(%s).*%e)./(%s)./%e./s.^2',Ho,lamda,phi,phi,r,phi,rw,gamma_2,Co);

[Walter Roberson]

why are you constructing string of those instead of symbolic expression ?

[Nat A]

To be honest I don't really know. I was told to use sprintf as above because it's easier to work with when using complicated equations. Most everything I'm using has been inherited from older things and adapted as needed by previous people and now myself and I don't doubt things have gotten mangled along the way. I'm not really fimiliar with the finer points of Matlab code so any help you could give to make this work better would be amazing.

[Nat A]

better?

syms(s)
phi = S_est./T_est.*s;
E_3 = (Co.*s.*besselk(0,sqrt(phi).*rw)+lamda.*sqrt(phi).*besselk(1,sqrt(phi).*rw))./(Co.*s);
E_4 = (Cs.*s.*besselk(0,sqrt(phi).*rw)+lamda.*sqrt(phi).*besselk(1,sqrt(phi).*rw))./(Cs.*s);
gamma_2 = besselk(0,sqrt(phi).*r).^2-(E_3).*(E_4);
h_ob_b = -Ho.*lamda.*sqrt(phi).*besselk(0,sqrt(phi).*r).*besselk(1,sqrt(phi).*rw)./(gamma_2)./Co./s.^2;

I now get the error:

Error using str2sym (line 46)
Argument must be string, character vector, or cell array of character vectors.

Error in INVLAP_2018_d (line 31)
 btf = matlabFunction(bt*str2sym(Fs), 'vars', {s, bt})

[Walter Roberson]

After that change you would just use that return value without str2sym as it is already sym. You would still use matlabFunction though.

[Nat A]

OK, thanks. What I now have is:

btf = matlabFunction((bt.*h_ob_b));
for kt=1:nnt                                
   tt=radt(kt);
   s=alfa/tt;                 
   bt=beta/tt; 
   test=double(subs(btf))
   ft(kt)=sum(real(test)) 
end

'test' converts the symbolic equation to a complex number, I can see this outputted when I run the code, but 'ft' is reported as a symbolic equation rather than a number. All the other variables within the 'for/end' loop are numeric. When the run is finished and I enter sum(real(test)) in the command window though a real number is returned. Do you know what the problem might be?

[Walter Roberson]

The implication would seem to be that somehow it was thinking that ft was initialized to symbolic before that point. I looked back over your code and code not see it, but perhaps I was not looking at the current version ?

[Nat A]

The code's been changing as we're working through it. ft was not previously defined but I did so explicitly now and it looks like it fixed the problem. It's still running too slow to be of much use however; previously execution with with nnt=10000 was ~4seconds. The current version is:

r=8; rw=0.0762; rcs = 0.0254; rco = 0.0254; Ho=2; S_est = 2e-6; T_est = 2e-5; lamda=9.5756e-06;
lamda = 2*pi()*rw*T_est; Cs = pi()*rcs^2; Co = pi()*rco^2;
a=6; ns=20; nd=19; tini=1; tend=10000; nnt=1000; radt=linspace(tini,tend,nnt); 

syms s
phi=S_est./T_est.*s;
E_3=(Co.*s.*besselk(0,sqrt(phi).*rw)+lamda.*sqrt(phi).*besselk(1,sqrt(phi).*rw))./(Co.*s);
E_4=(Cs.*s.*besselk(0,sqrt(phi).*rw)+lamda.*sqrt(phi).*besselk(1,sqrt(phi).*rw))./(Cs.*s);
gamma_2=besselk(0,sqrt(phi).*r).^2-(E_3).*(E_4);
h_ob_b=-Ho.*lamda.*sqrt(phi).*besselk(0,sqrt(phi).*r).*besselk(1,sqrt(phi).*rw)./(gamma_2)./Co./s.^2;

for n=1:ns+1+nd               
   alfa(n)=a+(n-1)*pi*j;
   beta(n)=-exp(a)*(-1)^n;
end

n=1:nd;
bdif=fliplr(cumsum(gamma(nd+1)./gamma(nd+2-n)./gamma(n)))./2^nd;
beta(ns+2:ns+1+nd)=beta(ns+2:ns+1+nd).*bdif;
beta(1)=beta(1)/2;
ft=zeros(nnt,1); syms bt %bt=0+0i;
btf = matlabFunction((bt.*h_ob_b));

for kt=1:nnt   
    kt
   tt=radt(kt);
   s=alfa/tt;                 
   bt=beta/tt; 
   test=double(subs(btf));
   ft(kt)=sum(real(test));
end

[Walter Roberson]

syms Alfa Beta TT
F = simplify(subs(bt.*h_ob_b, {s,bt}, {Alfa/TT, Beta/TT}));
FF = sum(real(subs(F,{Alfa,Beta},{alfa(:), beta(:)})));
FFF = matlabFunction(FF,'File','FFF.m','optimize',true);
ft = FFF(radt);

The creation of the .m file is not especially fast because of the optimization, but the execution for the 1000 points is less than 0.2 seconds.

You might also want to try timing with optimize turned off in the MATLAB function -- taking the tradeoff between the time spend optimizing the code and the time spend executing the optimized code.

[Nat A]

That's great, thanks so much!

Answer 1

[Walter Roberson]

(Reposting as an Answer)

syms Alfa Beta TT
F = simplify(subs(bt.*h_ob_b, {s,bt}, {Alfa/TT, Beta/TT}));
FF = sum(real(subs(F,{Alfa,Beta},{alfa(:), beta(:)})));
FFF = matlabFunction(FF,'File','FFF.m','optimize',true);
ft = FFF(radt);

The creation of the .m file is not especially fast because of the optimization, but the execution for the 1000 points is less than 0.2 seconds.

You might also want to try timing with optimize turned off in the MATLAB function -- taking the tradeoff between the time spend optimizing the code and the time spend executing the optimized code.