Calculation errors while using subs and det functions

9 次查看（过去 30 天）
Selçuk Sehitoglu
Selçuk Sehitoglu 2023-5-29
评论： Torsten 2023-5-31
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 !

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2023-5-30
编辑：Matt J 2023-5-30
Ran in:
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
  4 个评论
显示 2更早的评论
Selçuk Sehitoglu
Selçuk Sehitoglu 2023-5-31
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
Torsten 2023-5-31
Ran in:
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;

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by