Linear equations with symbolic matrix

I have a set of equations of the form Ax=b, but some of the elements in the matrix are also variables (one variable (y) which I have given as a symbolic variable, but is contained in multiple elements). I then have an additional condition (of the form ax1+bx2+cx3+dx4=e), which should make the system solvable. However, I tried adding this as an extra row to the matrix and then using x=linsolve(A,b), but this gives me an inconsistency error, presumably because I haven't told it I want to solve for the variable y. How can I go about combining these two things together?

Edit: I've tried making the x's symbolic variables as well, and then using vpasolve with eqn = A*x==b;, but I'm just getting empty syms as my output.

[Torsten]

If entries in the matrix A depend on a variable, you no longer have a linear system of equations and linsolve is not the adequate tool to solve. Use e.g "fsolve" instead.

[Alexander Holmes]

@Torsten Should the vpasolve method work?

[Torsten]

I can't understand why so many people prefer working with symbolic algebra if in the end, they want to get a numerical result. But yes - vpasolve should also work in your case.

[John D'Errico]

I'm also frequently flummoxed at the number of times I see people using symbolic variables for purely numerical problem. I think it has something to do with the design of MATLAB, in the sense that if you don't know the value of an unknown, then many people assume it must be a symbolic variable. Or it might be in how we teach people how to use MATLAB.

This is reflected in a simple example. Suppose I had a simple quadratic equation, f(x)=x^2-2. I'd like to solve for x such that it holds true. Assuming that I cannot do it using pencil and paper, in MATLAB, my choices would seem to be

1. Define x as a symbolic variable, then use solve.
2. Define a function of a variable x, then use a tool like fsolve or fzero.

Either is of course appropriate and totally valid.

In the first case, we do

syms x
xsol = solve(x^2 - 2,x)

xsol =
  2^(1/2)
 -2^(1/2)

In the second case, we might do it as

fun = @(x) x.^2 - 2;
xsol = fzero(fun,1)
xsol =
           1.4142135623731

But then I see people trying to use a tool like fzero or fsolve, but with symbolic objects. Of course, even though that can be done, it tends to be slower, and the symbolic variavbles were never neccesary at all.

Answer 1

[Sulaymon Eshkabilov]

Here are some of the ways of solving linear systems:

%% Given 13-by-13 system of linear equations
clearvars; clc
A = magic(13);                 % [A] Generated: with magic()
B = randi([-169,169], 13,1);   % [B] Generated: within [-169, 169] 

%% 1-Way: inv()  or pinv()     %% INVERSE matrix method 
x1a = inv(A)*B;
Err_INV = norm(A*x1a-B)/norm(B)     %#ok: ERROR checking 
x1b = pinv(A)*B;
Err_PINV = norm(A*x1b-B)/norm(B)     %#ok: ERROR checking 
%% 2-Way: \     %% backslash 
x2  = A\B;
Err_BACKSLASH = norm(A*x2-B)/norm(B)    %#ok: ERROR checking 
%% 3-Way: mldivide()          %% Left divide function 
x3 = mldivide(A, B);
Err_MLDIVIDE = norm(A*x3-B)/norm(B)   %#ok: ERROR checking 
%% 4-Way: Using linsolve()    
x4 = linsolve(A, B);
Err_LINSOLVE = norm(A*x4-B)/norm(B)   %#ok: ERROR checking 
%% 5-Way: Using lsqr()    
x5 = lsqr(A, B);
Err_LSQR = norm(A*x5-B)/norm(B)   %#ok: ERROR checking
%%   %%%%%%%%%%%%%%%%%%%%  Symbolic MATH   %%%%%%%%%%%%%%%%%%%%%%%%%%
%% Still there is a need to use symbolic MATH toolbox then
%%  solve()            %% SOLVE() symbolic math method 
x = sym('x', [1, 13]); x=x';
Eqn = A*(x); Eqn = Eqn - B;
Solution = solve(Eqn); SOLs = struct2array(Solution);
SOLs = double(SOLs); x13 = SOLs';
Err_SOLVE = norm(A*x13-B)/norm(B)   %#ok: ERROR checking

%%  VPASOLVE()
x = sym('x', [1, 13]); x=x';
Eqn = A*(x); Eqn = Eqn - B;
Solution = vpasolve(Eqn); SOLs = struct2array(Solution);
SOLs = double(SOLs); x13A = SOLs';
Err_VPASOLVE = norm(A*x13A-B)/norm(B) %#ok: ERROR checking

Note that among them VPASOLVE() is the most crude one and makes the largest error.

[Walter Roberson]

Note that mldivide is completely equivalent to \

The mathematics operators all have formal names that are invoked as functions, and the formal name of \ is mldivide.

For example, A+B*C invokes plus(A, mtimes(B, C)) . This is not an analogy, and not "just as if": when you invoke A+B then it literally converts it internally to plus(A,B) and does function resolution on the plus() function for the datatype of A and B and calls that.

[Walter Roberson]

vpasolve() does not the largest error. What has happened is that you have assumed the order of variables in the Solution returned by vpasolve(). vpasolve() happens to sort as x1, x10, x11, x12, x13, x2, x3, x4, x5, x6, x7, x8, x9 -- that is, it used sort() on the field names.

Answer 2

[Sulaymon Eshkabilov]

I don't see your point here. Why you are commenting all of these. In MATLAB, two types of commands that can be used interchangably, e.g.:

mldivide() vs. \ ,

mrdivide() vs. /,

ldivide() vs. .\

rdivide() vs. /,

transpose() vs. '

mtime() vs. *

...

etc.

[Walter Roberson]

Your code calculates Err_BACKSLASH and Err_MLDIVIDE separately, which could leave some people with the impression that the \ operator and mldivide() were different calculation methods, potentially leaving people confused about the difference between them. I added a note of clarification for readers that they are exactly the same, so as to reduce the potential for confusion.

[madhan ravi]

Second sir Walter’s comment Vs is not the same meaning as as .

Answer 3

[Sulaymon Eshkabilov]

Walter: your comments on vpasolve() are not relevant. pl., be kind and Run the code line by line and will understand how it computes the solutions and the computational error.

[Walter Roberson]

>> Solution = vpasolve(Eqn), SOLs = struct2array(Solution)

Solution =

struct with fields:

x1: [1×1 sym]

x10: [1×1 sym]

x11: [1×1 sym]

x12: [1×1 sym]

x13: [1×1 sym]

x2: [1×1 sym]

x3: [1×1 sym]

x4: [1×1 sym]

x5: [1×1 sym]

x6: [1×1 sym]

x7: [1×1 sym]

x8: [1×1 sym]

x9: [1×1 sym]

SOLs =

[ -1.6839000402052454032513604445969, -1.6606437889343776094903860519936, 1.3760517410961011558880708225448, -0.015793185219017528178139988945989, 0.54030441526167947929170470768636, 0.95168286561436454467679382384612, 0.10708655577622194915605654268839, -0.89983202171050798258438522615755, -0.028811615623429717726959907020098, -0.10828923358180420209389990696358, 1.6468300929839391377852916314455, 0.039134777785474043308821081081873, -0.14074364016647478985853016053839]

>> Solution.x10

ans =

-1.6606437889343776094903860519936

Notice that the second element of the struct2array() result, SOLs, is taken from the second field listed in SOLs, the name with field named x10. That is what struct2array() does: converts the entries by field order. struct2array() calls struct2cell() and then [] the cell entries together.

You then substitute that SOLs in for x values by position, which implies that you believe that SOLs(2) corresponds to x2, which it does not in this case.

>> Solution = solve(Eqn), SOLs = struct2array(Solution)

Solution =

struct with fields:

x1: [1×1 sym]

x2: [1×1 sym]

x3: [1×1 sym]

x4: [1×1 sym]

x5: [1×1 sym]

x6: [1×1 sym]

x7: [1×1 sym]

x8: [1×1 sym]

x9: [1×1 sym]

x10: [1×1 sym]

x11: [1×1 sym]

x12: [1×1 sym]

x13: [1×1 sym]

Notice that solve() does happen to return fields in the order x1, x2, x3, and so on, so there is no problem doing the struct2array() on that result.

When you take into account correct ordering of the fields you will see

>> Solution = vpasolve(Eqn); norm(subs(A*x,Solution)-B)/norm(B)

ans =

4.544317137566036944739238644881e-40

which is the smallest error other than the 0 error you can get with solve(). Furthermore, if you convert to double for fair comparison with your other possibilities, then

>> Solution = vpasolve(Eqn); norm(double(subs(A*x,Solution))-B)/norm(B)

ans =

0

which is as good as solve().

The difference is just down to correct ordering of the results.

Or perhaps you would prefer to see:

Solution = vpasolve(Eqn); SOLs = struct2array(Solution); SOLs = SOLs([1 6:13 2:5]);
SOLs = double(SOLs); x13A = SOLs';
Err_VPASOLVE = norm(A*x13A-B)/norm(B) %#ok: ERROR checking

Err_VPASOLVE =

2.12678054996351e-16

which is lower even than Err_SOLVE which is 2.57365786649384e-16