null(Matrix) returns an empty vector

54 次查看（过去 30 天）
Albert Heinrichs
Albert Heinrichs 2021-1-3
评论： John D'Errico 2021-1-3
I can't get the null space of a Matrix. This is how the matrix is calculated (use values Ta_in=3000, k_in=3, inputnumber=1 as reference):
function Mat = Matrix(Ta_in,k_in,inputnumber)
syms zeta(x) Psi(x) V(x) ypsilon0(x) sigma k Ta ypsilon0 Dypsilon0 mu
Ta=Ta_in;
k=k_in;
number=inputnumber;
mu_vec=[0,0,0,0,0.5,1];
mu=mu_vec(number);
ypsilon0(x)=1/(2*sigma+sigma^2)*((-1-sigma*x)+(1+sigma)^2/(1+sigma*x)+mu*(1+sigma)^2*(1+sigma*x-1/(1+sigma*x)));
ypsilon0_0=limit(ypsilon0(x),sigma,0);
Dypsilon0_0=limit(diff(ypsilon0(x),x),sigma,0);
sigma_vec=[0,0.5,1,2,0,0];
sigma=sigma_vec(number);
if sigma==0
ypsilon0=ypsilon0_0;
Dypsilon0=Dypsilon0_0;
else
ypsilon0=ypsilon0(x);
Dypsilon0=diff(ypsilon0,x,1);
end
ypsilon0=subs(ypsilon0);
Dypsilon0=subs(Dypsilon0);
gamma=0;
eqs=[zeta(x)==1/(1+sigma*x)*(diff(Psi(x),x,2)-k^2*Psi(x))-sigma*diff(Psi(x),x,1)/(1+sigma*x)^2,...
gamma*zeta(x)-Ta*2*ypsilon0*k*V(x)/(1+sigma*x)==diff(zeta(x),x,2)+sigma*diff(zeta(x),x,1)/(1+sigma*x)-(k^2+sigma^2/(1+sigma*x)^2)*zeta(x),...
gamma*V(x)-k*Psi(x)*Dypsilon0+sigma*ypsilon0/(1+sigma*x)/(1+sigma*x)==diff(V(x),x,2)+sigma/(1+sigma*x)*diff(V(x),x,1)-(k^2+sigma^2/(1+sigma*x)^2)*V(x)];
vars=[zeta(x), Psi(x), V(x)];
[newEqs,newVars]=reduceDifferentialOrder(eqs,vars);
newEqs=subs(newEqs);
xystruct=struct('x',[],'y',[]);
f=zeros(3,3);
for n=1:3
switch n
case 1
c=[1 0 0];
case 2
c=[0 1 0];
case 3
c=[0 0 1];
end
%%initial conditions:
Psi0=0;
DPsi0=0;
zeta0=c(1);
Dzeta0=c(2);
V0=0;
DV0=c(3);
%%order: newVars =zeta(x), Psi(x), V(x), Dzeta_t(x), DPsit(x), DVt(x)
initConditions=[zeta0 Psi0 V0 Dzeta0 DPsi0 DV0];
[M,F]=massMatrixForm(newEqs,newVars);
fun = M\F;
odefun=odeFunction(fun,newVars);
opts = odeset('RelTol',1e-5,'AbsTol',1e-7);
[x,y]=ode45(odefun,[0 1], initConditions,opts);
xystruct(n).x=x;
xystruct(n).y=y;
[rows,cols]=size(y);
f(1,n)=y(rows,2);%%f1 = psi
f(2,n)=y(rows,5);%%f2 = Dpsi
f(3,n)=y(rows,3);%%f3 = V
end
Mat=f;
Then:
null(Matrix(3000,3,1))
This returns an empty vector.
Any ideas how to solve this issue?

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

采纳的回答

John D'Errico
John D'Errico 2021-1-3
编辑：John D'Errico 2021-1-3
That null returns an empty vector merely means your matrix is full rank. You cannot compute a non-empty null space of a full rank matrix. It does not exist.
A good question is to compute the SVD of your matrix. you need only the singular values. Is one of them close to being zero, but just a bit more than eps times the largest singular vlaue?
Mat
Mat =
0.724126524459701 0.285725197854218 -173.130875191614
-0.488098140173286 0.942936238526739 -938.432889302303
0.244437311796916 0.0428476702924243 -9.40505165295999
>> rank(Mat)
ans =
3
>> svd(Mat)
ans =
954.31657322839
0.845977933126111
0.000784006154677659
So not even that close to being singular. It has one moderately small singular value. If, in fact, your code is correct, but there are numerical issues corrupting the problem, the closest thing this matrix has to a null space vector is V(:,3) below:
[U,S,V] = svd(Mat);
V(:,3)
ans =
0.135712466140226
-0.990747692346457
-0.00106613522064244
Is it truly a nullspace vector? Well, no. Not very close.
Mat*V(:,3)
ans =
-0.000229246749867018
3.47877296913435e-05
0.000748932070224632
but as close as you can get.
  2 个评论
Albert Heinrichs
Albert Heinrichs 2021-1-3
ok thank you very much for your quick response! Would it help having the determinant of the Matrix be very close to zero (using svd)? because that's the point where i want to calculate/approximate the null space vector anyways. I've already iteratively approximated several points of interest e.g.
Ta=1.696240234375000e+03
k=3.2
n=1
John D'Errico
John D'Errico 2021-1-3
Ran in:
A matrix with determinant very close to zero. Sorry, but regardless what you learned from some class, the determinant is a poor way to learn anything about a matrix.
A = eye(10);
det(A)
ans = 1
det(10*A)
ans = 1.0000e+10
det(0.1*A)
ans = 1.0000e-10
Multiply or divide it by 10, and the determinant is either tiny or huge. But none of those matrices is even remotely singular.
Anyway, I already showed you how to find the closest thing to a null space vector from that matrix.
If you want to perturb the paramters to find a matrix that is singular, then instead of trying to push the determinant to zero, you will be better off trying to make the condition number of the matrix as large as possible.
help cond
Or, you can try to make the reciprocal condition number as small as possible.
rcond(Mat)
ans =
6.49202453580327e-07

请先登录，再进行评论。

更多回答（0 个）

类别

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

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by