Incorrect Substitution with subs() from Workspace Variables?

syms alpha s H t
system = 1/(s + 1/alpha);
y = ilaplace(1/s * (1-exp(-s*H)) * system)

y = 

Define specific values for alpha and H

% values
alpha = 15;
H = 32;

Substitute those values into y

foo = subs(y)

foo = 

Now substitute a specific value of t.

subs(foo,t,50)

ans = 

That looks like the correct result.

Define t as a workspace variable, and then try to subs in alpha, H, t all together, yields an incorrect result?

t = 50;
subs(y)

ans = 

0

Bug? Or am I misunderstanding something about @doc:subs?

Answer 1

[Piyush Kumar]

@Paul,

When you substitute all variables at once in MATLAB using subs(y, [alpha, H, t], [15, 32, 50]), MATLAB might not handle the substitutions in the expected order, leading to an incorrect result.

However, when you substitute variables one by one, MATLAB updates the expression incrementally, ensuring that each substitution is correctly applied in sequence. This avoids any conflicts or misinterpretations that might occur when substituting multiple variables at once.

Ensuring Correct Substitution

To ensure correct substitution, you can continue using the step-by-step approach, which helps MATLAB correctly map each variable:

% Define symbolic variables
syms alpha s H t

% Define the system
system = 1/(s + 1/alpha);
y = ilaplace(1/s * (1-exp(-s*H)) * system);

% Substitute specific values for alpha and H
alpha_value = 15;
H_value = 32;
t_value = 50;

% Substitute values one by one
foo = subs(y, alpha, alpha_value);
foo = subs(foo, H, H_value);
foo = subs(foo, t, t_value);

% Display the result
disp('Result after step-by-step substitution:');

Result after step-by-step substitution:

disp(foo);

Let's consider an example where the order of substitution can lead to different results.

Consider the expression:

Now, let's substitute x = y + z, y = z + 1, and z = 2.

Substitution Order 1:

Substitute z = 2 =>

Substitute y = z + 1 =>

Substitute x = y + z =>

In this case, the expression becomes undefined due to division by zero.

Substitution Order 2:

Substitute x = y + z =>

Substitute y = z + 1 =>

Substitute z = 2 =>

In this case, the expression evaluates to 5.

The order of substitution can significantly affect the final result, especially in complex expressions where variables are interdependent. In the first order, the expression becomes undefined due to division by zero, while in the second order, the expression evaluates to a finite value.

​

[Paul]

At first I wasn't sure how this answer applies to the problem, but it got me thinking about ordering of the substitutions.

Define the expression.

syms alpha s H t
system = 1/(s + 1/alpha);
y = ilaplace(1/s * (1-exp(-s*H)) * system)

y = 

Now do the subs of alpha, H, and t in all six possible orders.

y1 = subs(y,alpha,15);y1 = subs(y1,H,32); y1 = subs(y1,t,50)

y1 = 

y1 = subs(y,alpha,15);y1 = subs(y1,t,50); y1 = subs(y1,H,32)

y1 = 

0

y1 = subs(y,H,32);y1 = subs(y1,alpha,15); y1 = subs(y1,t,50)

y1 = 

y1 = subs(y,H,32); y1 = subs(y1,t,50);y1 = subs(y1,alpha,15)

y1 = 

y1 = subs(y,t,50); y1 = subs(y1,alpha,15);y1 = subs(y1,H,32)

y1 = 

0

y1 = subs(y,t,50); y1 = subs(y1,H,32);y1 = subs(y1,alpha,15)

y1 = 

0

We see that the incorrect result obtains when t is subbed before H.

Define the ilaplace part of the expression.

Y = ilaplace(exp(-H*s)/s/(s+1/alpha),s,t)

Y = 

If we subs H first we get an expression in t

subs(Y,H,32)

ans = 

●

And then we sub t to get the final result

subs(ans,t,50)

ans = 

If we subs t first, we get an ilaplace expression with '50' as the transformation variable.

subs(Y,t,50)

ans = 

At first glance, that looks a bit odd. But, based on my experience with fourier, I'd assume that expression means to find the ilaplace of the first argument and then sub in the value of 50 for transformation variable.

But, when we sub in the value of H we get

subs(ans,H,32)

ans = 

which is what we'd get for the ilaplace if H == 0

ilaplace(1/s/(s+1/alpha),s,50)

ans = 

So now it looks like a bug in ilaplace.

Note that even if we define the exponential term explicitly, we still get the wrong answer.

ilaplace(exp(-32*s)*1/s/(s+1/alpha),s,50)

ans = 

whereas the correct answer is

subs(ilaplace(exp(-32*s)*1/s/(s+1/alpha),s,t),t,50)

ans = 

So still looks like a bug in ilaplace. If one is allowed to use a constant as the third argument, ilaplace should return the correct result. If that's not supposed to be allowable, then it should throw an error.

[Paul]

MathWorks Tech Support said that they will consider fixing ilaplace in a future release.

[Walter Roberson]

When you substitute all variables at once in MATLAB using subs(y, [alpha, H, t], [15, 32, 50]), MATLAB might not handle the substitutions in the expected order, leading to an incorrect result.

When you subs() using a vector of variable names and a vector of values, subs() reads off the variable names and the corresponding replacement values, and does the substitutions effectively simultaneously. For example

syms a b
f = 10*a + b

f = 

subs(f, [a b], [b a])

ans = 

This contrasts with iterative order, which is not used

subs(subs(f, a, b), b, a)

ans = 

subs(subs(f, b, a), a, b)

ans = 

Simultaneous order is also necessary to get matrix subsitution right

subs([a, b], {a, b}, {[1;2], [3;4]})

ans = 

subs(subs([a, b], {a}, {[1;2]}), {b}, {[3;4]})

ans = 

subs(subs([a, b], {b}, {[3;4]}), {a}, {[1;2]})

ans = 

Therefore, doing simutaneous substitution does produce the expected results -- unless, that is, what is expected is wrong.

subs(10*a + b, [a b], [b 3])

ans = 

If the expectation was that the 10*a would first be processed through a->b, giving 10*b+b --> 11*b, and that then the b->3 would be processed, giving 33 as the final result, then Yes, the actual results might be "unexpected"

But when something like

subs(y, [alpha, H, t], [15, 32, 50])

is processed, the order of simultaneous substitutions is irrelevant.

The fact that the substitutions is not iterative can in theory matter for situations such as diff(f,alpha) as iterative substitution could end up turning the diff(f,alpha) into diff(f,15) which would be a request for the 15th derivative of the default variable.