[問題]3000P~~為什麼矩陣明明有值 算出來卻是NaN?已刪文

看板MATLAB作者 (綠扁帽)時間6年前 (2018/05/17 16:50), 6年前編輯推噓6(606)
留言12則, 3人參與, 6年前最新討論串1/1
% Kalman filter for the example of free body falling % Using previous data to modified clear all clc I=eye(360,360); N=160; y=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 4 6 9.8 10 10.2 17.8 26 29 33.4 35.6 53.4 61 70 75.6 81.1 82 82.2 79.8 77.8 79 80.8 79 77.8 71.5 65 62.2 58 55 50.5 48.7 47.7 45.5 42.2 37.8 37.8 34.5 33.4 28.8 25.5 23.3 22.2 21.1 20 17.8 17.1 16.4 15.7 15.5 14.5 13.4 12.2 11.2 10 9.5 9 9.5 10 8.9 5.4 5.5 5.6 5.5 5.4 5.3 5.2 5.1 4.4 4 4 4 4 4 4 3 2 1 0 0.7 1.4 2.2 1.4 0.7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] dt=60 gg=[19.26 -10.818 0]; hh=eye(1,360); aa=conv2(hh,gg,'same') gg=[-6.969 11.818 -3.849]; hh=eye(358,360); bb=conv2(hh,gg,'same') cc=[0 0 0 0 0 . . . . . 0 -3.12 4.12] % cc大小為(1,360) C=[aa;bb;cc] %將aa、bb、cc組合成C矩陣 U=0 M = eye(360) A = (M/dt)^-1*((M/dt)-C); U = 0 ; H = zeros(1,360); H(1,1) = 1; x=zeros(360,1,N); x(:,1,1)=[575831 0 0 0 0 . . . . 0]; % x(:,1,1)是360*1的矩陣 x0=zeros(1,N) p=zeros(360,360,N); p(:,:,1)=eye(360); F=zeros(360,1); for k=2:N x(:,1,k) = A*x(:,:,k-1)+U; p(:,:,k) = A*p(:,:,k-1)*A'; F(:,1,k) = p(:,:,k)*H' / (H*p(:,:,k)*H'); x(:,1,k)=x(:,:,k)+F(:,1,k)*(y(k-1)-H*x(:,:,k-1)); x0(1,k)=x(230,1,k); p(:,:,k)=(I-F(:,1,k)*H)*p(:,:,k);; end y x0 t=1:N; plot(t,y,t,x0); ===================================== 我要求的是x0矩陣的結果,但算出來卻是很多NaN... 但是矩陣用手算的話應該是有值的... 昨天查一下 好像有的時候算矩陣會出現這種結果 (可能是分母太接近0的關係??) 也可能是矩陣算出來太大或太小的原因? 不知道要不要換成sym去算=.=.. 跪求高手教學... -- ※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 125.230.66.216 ※ 文章網址: https://www.ptt.cc/bbs/MATLAB/M.1526547042.A.48F.html ※ 編輯: GreenBeret (125.230.66.216), 05/17/2018 16:52:27

05/17 17:10, 6年前 , 1F
先看在哪個k出問題,我猜是F = ... / ... 除到零了
05/17 17:10, 1F

05/17 20:16, 6年前 , 2F
singular matrix 檢查eigenvalues
05/17 20:16, 2F

05/18 17:37, 6年前 , 3F
x(:,1,k) = A*x(:,:,k-1)+U; 就出問題了... 一堆NaN
05/18 17:37, 3F

05/18 19:26, 6年前 , 4F
那就看一下第一次算A的時候對不對
05/18 19:26, 4F

05/18 23:43, 6年前 , 5F
感謝 我試看看
05/18 23:43, 5F

05/19 16:50, 6年前 , 6F
等等 A = (M/dt)^-1*((M/dt)-C); ^-1不是反矩陣啊
05/19 16:50, 6F

05/19 16:53, 6年前 , 7F
其實沒差 M是diagonal
05/19 16:53, 7F

05/19 16:55, 6年前 , 8F
阿 不對 這樣0就NaN了 把^-1改成inv()應該就可以
05/19 16:55, 8F

05/20 01:42, 6年前 , 9F
結果inv()還是沒辦法XD 我再查一下到底出錯在哪裡
05/20 01:42, 9F

05/20 01:43, 6年前 , 10F
3Q
05/20 01:43, 10F

05/20 01:50, 6年前 , 11F
我早上算第一次的A的時候....覺得是因為太接近0了Q_Q
05/20 01:50, 11F

05/20 01:50, 6年前 , 12F
我想個辦法
05/20 01:50, 12F
文章代碼(AID): #1Q_K9YIF (MATLAB)
文章代碼(AID): #1Q_K9YIF (MATLAB)