How to Solve a System of Equations for symfun Objects?

Suppose I have a system of equations:

syms R E Y G
eq1 = Y == G*E;
eq2 = E == R-Y;

Solve for Y and E in terms of R

sol = solve([eq1,eq2],[Y,E]);

Make an equation to show the solutions for Y and E

eqn = [Y == sol.Y;E == sol.E]

eqn = 

But Y, R, G, and E are all really functions of s, so

syms Y(s) R(s) G(s) E(s)
eqn = subs(eqn)

eqn = 

1. Why did E in eqn(2) not get converted to E(s)?

2. Is there a way to define all of the R(s),E(s),G(s),Y(s) as symfun from the start and simultaneously solve eq1 and eq2 directly for Y(s) and E(s)? @doc:solve doesn't work (or at least I couldn't figure it out) with equations of symfun objects, and @doc:isolate only works with scalar expressions. I looked at using @doc:mapSymType to temporarily convert the symfun objects in eq1 and eq2 to sym, then @doc:solve, then convert back to symfun, but that path looked a bit tortuous.

[Torsten]

1. Replacing E by F works :-)

The display style for E also differs from Y, G and R. Maybe a reserved constant or function name ? But if it existed, it should be overwritten by the "syms" definition.

[Paul]

I thought something looked strange about that equation, but did not notice the E was displayed differently. But, to be clear, it's not just a display thing.

syms R E Y G
eq1 = Y == G*E;
eq2 = E == R-Y;
sol = solve([eq1,eq2],[Y,E]);
eqn = [Y == sol.Y;E == sol.E];
syms Y(s) R(s) G(s) E(s)
eqn = subs(eqn)

eqn = 

symType(lhs(eqn))

ans = 2x1 string array
    "symfun"
    "variable"

Apparently the subs function works differently for E/E(s). As you said, maybe it is a reserved object of some sort. Seems like a bug to me, vaguely (but not really) similar to this thread.

[Walter Roberson]

syms A B C D E F G H I J K L M N O P Q R S T U V W X Y Z
syms a b c d e f g h i j k l m n o p q r s t u v w x y z
[A B C D E F G H I J K L M N O P Q R S T U V W X Y Z]

ans = 

[a b c d e f g h i j k l m n o p q r s t u v w x y z]

ans = 

It looks like there might be special treatment for D E I O

eq1 = D == G*E;
eq2 = E == R-Y;
eq3 = I == E+Y;
eq4 = O == R*Y;
eq5 = N == G+E;

sol = solve([eq1,eq2,eq3,eq4,eq5],[D,E,I,O,N])

sol = struct with fields:
    D: G*(R - Y)
    E: R - Y
    I: R
    O: R*Y
    N: G + R - Y

eqn = [D == sol.E;E == sol.E;I == sol.I;O == sol.O;N == sol.N]

eqn = 

syms D(s) E(s) I(s) N(s) O(s) G(s) R(s)
subs(eqn)

ans = 

●

Observe that all of D E I O do not substitute properly, but N (not on the suspect list) does substitute properly

[Torsten]

Strange that no one has noticed this yet.

[Paul]

syms A B C D E F G H I J K L M N O P Q R S T U V W X Y Z
syms a b c d e f g h i j k l m n o p q r s t u v w x y z
eq1 = D == G*E;
eq2 = E == R-Y;
eq3 = I == E+Y;
eq4 = O == R*Y;
eq5 = N == G+E;
sol = solve([eq1,eq2,eq3,eq4,eq5],[D,E,I,O,N])

sol = struct with fields:
    D: G*(R - Y)
    E: R - Y
    I: R
    O: R*Y
    N: G + R - Y

eqn = [D == sol.E;E == sol.E;I == sol.I;O == sol.O;N == sol.N]

eqn = 

D = d;
E = e;
I = i;
O = o;
subs(eqn)

ans = 

●

subs seems to properly substitute garden variety sym objects into D/E/I/O, just not symfun objects. Still seems like bug to me (display and subs).

Also, so as to not lose the forest for the trees, I'm more interested in item (2) in the original question.

[Paul]

I submitted a bug report. Response is that developers will consider a fix for a future Matlab release, presumably for both the display and subs issue for D, E, I, and O

Answer 1

[Walter Roberson]

Is there a way to define all of the R(s),E(s),G(s),Y(s) as symfun from the start and simultaneously solve eq1 and eq2 directly for Y(s) and E(s)?

No, there is not.

[Paul]

That's too bad. I couldn't even find a way to convert an unassigned symfun to a sym, at least not with reasonable effort.

For this particular problem, I'm pretty sure there is a path by converting eq1 and eq2 to strings using char(), replacing the (s) with "", then convert back with str2sym, then solve(), the subs to get it all back to symfuns, with some other intermediate steps needed as well.

[Walter Roberson]

syms R(s) e(s) Y(s) G(s) H(s) F(s)  % use lower case 'e' to avoid problems discussed above with 'E'
eqin(1,1) = Y(s) == G(s)*e(s);
eqin(2,1) = e(s) == R(s) - F(s);
eqin(3,1) = F(s) == H(s)*Y(s);
eqsf = findSymType(eqin, 'symfun');
eqsyms = arrayfun(@(SF) feval(symengine, 'op', SF, 0), eqsf);
eqsubs = subs(eqin, eqsf, eqsyms)

eqsubs = 

lhssyms = arrayfun(@(SF) feval(symengine, 'op', SF, 0), findSymType(lhs(eqin), 'symfun'));
subs(lhssyms == subs(lhssyms, solve(eqsubs, lhssyms)), eqsyms, eqsf)

ans = 

Answer 2

[Paul]

I took a shot. Here's a start, doesn't do any error checking, assumes that there is one valid solution ....

syms R(s) e(s) Y(s) G(s) H(s) F(s)  % use lower case 'e' to avoid problems discussed above with 'E'
eqin(1,1) = Y(s) == G(s)*e(s);
eqin(2,1) = e(s) == R(s) - F(s);
eqin(3,1) = F(s) == H(s)*Y(s);

eqin

eqin = 

eqout = solvesymfun(eqin,{Y,e,F})

eqout = 

function eqout = solvesymfun(eqin,solfuns)

% argument of the symfuns in the equations to be solved
varg = argnames(solfuns{1});

% recast the equation in terms of sym rather than symfun
eqout = str2sym(replace(char(eqin),"("+char(varg)+")",""));

% variables to solve for, as both string and sym
for ii = 1:numel(solfuns)
    v2str(ii) = string(replace(char(solfuns{ii}),"("+char(varg)+")",""));
    v2sym(ii) = sym(v2str(ii));
end

% solve
sol = solve(eqout,v2sym);

% create output equation
for ii = 1:numel(solfuns)
    eqout(ii) = v2sym(ii) == sol.(v2str(ii));
end

% return to symfun
eqout = subs(eqout);

end