Calculation errors while using subs and det functions

Hello everybody, i have a 9*9 symbolic matrix and i need to calculate thatof determinant in an interval. I would like to insert a numeric array into symbolic variable and get the determinant values and i want to plot them at the end.

The problems are;

1- when i insert numerical array using "subs" function, then get the numerical expression using the "double" function, then calculate the determinant using the "det" function, there are big calculation errors.

syms beta_sym

beta_ara = [(0.0001:0.0001:0.001)';(0.002:0.001:0.01)';(0.02:0.01:1)';(1:1:20)';(25:5:200)'];
f_beta_num=zeros(length(beta_ara),1);
ddel = besselj(1,beta_sym)*beta_sym^8; %% "ddel" is an expression including
%% bessel & mod. bessel functions and 8th
%% order polynomial of beta_sym and so on

for i = 1:length(beta_ara)
    f_beta_num(i) = det(double(subs(ddel,beta_sym,beta_ara(i))));
end

figure(1);clf;
plot(beta_ara,(f_beta_num),'LineWidth',2);hold on;

2- when i use "det" function at the most inner part of the row it takes very long time, thus i don't prefer this alternative.

There must be something missing at the step where it calculates determinant after inserting the numerical value. I guess there must be a row added in order to correct calculation / calculation method.

Following you can see how the plot supposed to be & how it is:

Thanks in advance !

Answer 1

[Matt J]

You haven't given us access to the input variables needed to repeat the computation. However, the determinant you are evaluating is an order 9 polynomial function of the matrix entries. An order 9 polynomial will have very steep sections where small numerical errors make a big difference in the determinant (as illustrated below). It is usually a bad idea to use determinants for any numerical work. Ordinarily, you would use rcond.

A=rand(9)*diag(linspace(0,200,9))*rand(9);

det(A)

ans = 1.9458

det(A+rand(9)/1e9) %A with small errors

ans = 1.1326e+04

[Selçuk Sehitoglu]

Thanks for your answer. I have added the first row of the code which was absent and uncomment the row which includes "ddel". Unfortunately i have to plot the graph which i have added to question using the numerical value of the determinant. And i also should mention that the following code works perfect but takes too much time since i will use for-loop.

syms beta_sym
beta_ara = [(0.0001:0.0001:0.001)';(0.002:0.001:0.01)';(0.02:0.01:1)';(1:1:20)';(25:5:200)'];
f_beta_num=zeros(length(beta_ara),1);
ddel = besselj(1,beta_sym)*beta_sym^8; %% "ddel" is an expression including
%% bessel & mod. bessel functions and 8th
%% order polynomial of beta_sym and so on
for i = 1:length(beta_ara)
    f_beta_num(i) = double(subs(det(ddel),beta_sym,beta_ara(i)));
end
figure(1);clf;
plot(beta_ara,(f_beta_num),'LineWidth',2);hold on;

[Torsten]

I don't understand why you use "det" in the above code.

ddel is a scalar value - thus det(ddel) = ddel.

You will get the same results if you use

f_beta_num(i) = double(subs(ddel,beta_sym,beta_ara(i)));

instead of

f_beta_num(i) = double(subs(det(ddel),beta_sym,beta_ara(i)));

[Selçuk Sehitoglu]

I have maken a mistake when i rewrite the code, corrected form is:

syms beta_sym c11 c12
beta_ara = [(0.0001:0.0001:0.001)';(0.002:0.001:0.01)';(0.02:0.01:1)';(1:1:20)';(25:5:200)'];
f_beta_num=zeros(length(beta_ara),1);
denk1 = c11*besselj(1,beta_sym) + c12*besselk(1,beta_sym) + beta_sym^8;
denk2 = c11*beta_sym*besseli(1,beta_sym) + c12*bessely(1,beta_sym)*beta_sym + beta_sym^4;
denk = [denk1; denk2];
c = [c11; c12];
ddel = equationsToMatrix(denk, c);
for i = 1:length(beta_ara)
    f_beta_num(i) = double(subs(det(ddel),beta_sym,beta_ara(i)));
end
figure(1);clf;
plot(beta_ara,(f_beta_num),'LineWidth',2);hold on;

which works fine but slow.

And here is the code:

syms beta_sym c11 c12
beta_ara = [(0.0001:0.0001:0.001)';(0.002:0.001:0.01)';(0.02:0.01:1)';(1:1:20)';(25:5:200)'];
f_beta_num=zeros(length(beta_ara),1);
denk1 = c11*besselj(1,beta_sym) + c12*besselk(1,beta_sym) + beta_sym^8;
denk2 = c11*beta_sym*besseli(1,beta_sym) + c12*bessely(1,beta_sym)*beta_sym + beta_sym^4;
denk = [denk1; denk2];
c = [c11; c12];
ddel = equationsToMatrix(denk , c);
for i = 1:length(beta_ara)
    f_beta_num(i) = det(double(subs(ddel,beta_sym,beta_ara(i))));
end
figure(1);clf;
plot(beta_ara,(f_beta_num),'LineWidth',2);hold on;

which works fast but problematic.

Sorry for the mistakes...

[Torsten]

syms beta_sym c11 c12
beta_ara = [(0.0001:0.0001:0.001)';(0.002:0.001:0.01)';(0.02:0.01:1)';(1:1:20)';(25:5:200)'];
f_beta_num=zeros(length(beta_ara),1);
denk1 = c11*besselj(1,beta_sym) + c12*besselk(1,beta_sym) + beta_sym^8;
denk2 = c11*beta_sym*besseli(1,beta_sym) + c12*bessely(1,beta_sym)*beta_sym + beta_sym^4;
denk = [denk1; denk2];
c = [c11; c12];
ddel = equationsToMatrix(denk, c);
det_ddel = matlabFunction(det(ddel));
f_beta_num = det_ddel(beta_ara);
figure(1);clf;
plot(beta_ara,f_beta_num,'LineWidth',2);hold on;