Re: [問題] polyfit已回收

看板MATLAB作者hilosi時間17年前 (2008/11/22 21:26)推噓1(1推 0噓 0→)

留言1則, 1人參與討論串3/5 (看更多)

※ 引述《AllenTing (阿丁)》之銘言： : ※ 引述《AllenTing (阿丁)》之銘言： : : 處理實驗數據的時候 : : 圖形經常都是y = ax^(-0.XXXX) 等的趨勢(EXCEL求出的趨勢線) : : 用matlab的Polyfit 都只能逼出-1,-2,-3次的方程式 : : 導致之後要求此圖形的斜率圖會很不準 : : 不知道有沒有其他方法可以順利的逼出精準的方程式 : : 不然實驗數據好難處理喔= = : : 謝謝 : 這是我要分析的數據 : d0 = 0.9956 : d = 0.9953 0.9953 0.9950 0.9948 0.9941 0.9906 0.9860 : r = 69.8333 66.6667 62.3333 59.4333 52.9000 36.0667 26.9000 : m = 0.0050 0.0100 0.0250 0.0500 0.1000 0.4000 0.8000 : 要做 m-r的圖 : 然後求出每一點 : Ma = 74.2; : whi = (Ma - 1000./m.*(d-d0)./d0)./d; : 求 Va = whi + m.*(d(whi)/d(√m)); : u = -1/8.314/303.15.*((d(r)/d(√m)); : 以上兩個綠色的地方就是做不出來的地方 : 逼近m-whi的地方就逼不出來了更別說要求斜率下面一串程式碼給你參考可以求出 Cubic Spline 所需要使用的係數如果你有兩串資料為 X(1, :), Y(1, :) 那麼你使用那個 function 求出來的 CIcoes(2, :) 就是各個 X 點上所對應的 dY/dX 值的一次微分數值如果你的資料是 X(:, 1), Y(:, 1) 那麼你使用那個 function 求出來的 CIcoes(:, 2) 就是各個 X 點上所對應的 dY/dX 值的一次微分數值你想要求的東西應該可以用下面的指令呼叫下面的 function 求出 [ CIcoes ] = Cubic_Int ( 2或3或4或6, m資料數, sqrt(m), whi, 0, 0 ) dwhi_dsqrt(m)=CIcoes(2, :); 或 dwhi_dsqrt(m)=CIcoes(:, 2); 如果要內插求其他非 X 點上的 dY/dX 請你去找有關 Cubic Spline 的書看看該如何去使用這 function 求出的係數請將下面的程式碼存在檔名為 Cubic_Int.m 的文字檔案中 -- function [ CIcoes ] = Cubic_Int ( CItype, CILen, CIX, CIY, CIm1, CImN ) %--------------------------------------------------------- % Cubic spline interpolation % Output : CIcoes % Input : CItype, Len, CIX, CIY, CIm1, CImN % % CIcoes : The Coefficients of Cubic Spline % the style of arrangement will be the same as X. % % % CItype : The types of Endpoint constraints % % 1. Clamped cubic spline, specify f'(X1), f'(Xn). % ( the best choice if the derivatives are known ) % % 2. Natural cubic spline. f''(X1)=0, f''(Xn)=0 % ( a relaxed curve ) % % 3. Extrapolate f''(X) to the endpoints. % % 4. f''(X) is constant near the endpoints. % f''(X1)=f''(X2), f''(Xn)=f''(Xn-1) % (Parabolically terminated Spline) % % 5. Specify S''(x) at each endpoint % % 6. Period Spline : f''(X1)=f''(Xn-1) % f''(Xn)=f''(X2) % % CILen : The total number of the discrete points. % % CIX : The x-magnitude of each discrete point are in order. % % CIY : The y-magnitude of each discrete point are in order. % % CIm1 : The endpoint constraint of the 1st discrete point. % % CImN : The endpoint constraint of the Nth discrete point. % % TDMA system equation matrices % % [b, c, 0, 0, 0] [f2] [d] % [a, b, c, 0, 0] [f2] [d] % ... ... = ... % [0, 0, a, b, c] [f2] [d] % [0, 0, 0, a, b] [f2] [d] %------- check the style of arrangement of CIX if length(CIX(1,:))==CILen CIX=CIX'; CIarr=1; else CIarr=0; end if length(CIY(1,:))==CILen CIY=CIY'; end CIf2=zeros(1, CILen); CIh=zeros(1, CILen-1); CID=zeros(1, CILen-1); CIu=zeros(1, CILen-2); for CIiter=1:CILen-1 CIh(1, CIiter)=CIX(CIiter+1, 1)-CIX(CIiter, 1); if CIh(1, CIiter)==0 fprintf('Find the same discrete point !!\n') else CID(1, CIiter)=(CIY(CIiter+1, 1)-... CIY(CIiter, 1))/CIh(1, CIiter); end end for CIiter=1:CILen-2 CIu(1, CIiter)=6*(CID(1, CIiter+1)-CID(1, CIiter)); end CIa=CIh(1, 2:CILen-2); CIb=2*(CIh(1, 1:CILen-2)+CIh(1, 2:CILen-1)); CIc=CIh(1, 2:CILen-2); CId=CIu(1, 1:CILen-2); CIup=zeros(1,CILen-4); switch CItype case 1 CIb(1,1)=CIb(1,1)-CIh(1,1)/2; CIb(1,CILen-2)=CIb(1,CILen-2)-... CIh(1,CILen-1)/2; CId(1,1)=CId(1,1)-3*(CID(1,1)-CIm1); CId(1,CILen-2)=CId(1,CILen-2)-... 3*(CImN-CID(1, CILen-1)); case 2 case 3 CIb(1,1)=CIb(1,1)+CIh(1,1)+... CIh(1,1)*CIh(1,1)/CIh(1,2); CIb(1,CILen-2)=CIb(1,CILen-2)+... CIh(1,CILen-1)+... CIh(1,CILen-1)*... CIh(1,CILen-1)... /CIh(1,CILen-2); CIc(1,1)=CIc(1,1)-CIh(1,1)*CIh(1,1)/CIh(1,2); CIa(1,CILen-3)=CIa(1,CILen-3)-... CIh(1,CILen-1)*... CIh(1,CILen-1)/... CIh(1,CILen-2); case 4 CIb(1,1)=CIb(1,1)+CIh(1,1); CIb(1,CILen-2)=CIb(1,CILen-2)+... CIh(1,CILen-1); case 5 CId(1,1)=CId(1,1)-CIh(1,1)*CIm1; CId(1,CILen-2)=CId(1,CILen-2)-... CIh(1,CILen-1)*CImN; case 6 CIup(1,1)=CIh(1,1); CIlow=CIh(1,CILen-1); end %---------------- Solving the TDMA ------------------- if CItype==6 CIb(1,CILen-2)=CIb(1,CILen-2)-... CIup(1,1)*CIlow/CIb(1,1); CId(1,CILen-2)=CId(1,CILen-2)-... CId(1,1)*CIlow/CIb(1,1); for CIiter=1:CILen-5 CIup(1,CIiter+1)=-CIa(1,CIiter)*... CIup(1,CIiter)/CIb(1,CIiter); end CIc(1,CILen-3)=CIc(1,CILen-3)-... CIa(1,CILen-4)*... CIup(1,CILen-4)/... CIb(1,CILen-4); for CIiter=1:CILen-3 CIb(1,CIiter+1)=CIb(1,CIiter+1)-... CIa(1,CIiter)*CIc(1,CIiter)/... CIb(1,CIiter); CId(1,CIiter+1)=CId(1,CIiter+1)-... CIa(1,CIiter)*... CId(1,CIiter)/... CIb(1,CIiter); end for CIiter=1:CILen-2 if CIiter==1 CIf2(1, CILen-CIiter)=... CId(1,CILen-CIiter-1)/... CIb(1,CILen-CIiter-1); elseif CIiter==2 CIf2(1, CILen-CIiter)=... (CId(1,CILen-CIiter-1)-... CIc(1,CILen-CIiter-1)*... CIf2(1, CILen-CIiter))/... CIb(1,CILen-CIiter-1); else CIf2(1, CILen-CIiter)=... (CId(1,CILen-CIiter-1)-... CIc(1,CILen-CIiter-1)*... CIf2(1, CILen-CIiter)-... CIup(1,CILen-CIiter-1)*... CIf2(1,CILen-1))/... CIb(1,CILen-CIiter-1); end end else for CIiter=1:CILen-3 CIb(1,CIiter+1)=CIb(1,CIiter+1)-... CIa(1,CIiter)*CIc(1,CIiter)/... CIb(1,CIiter); CId(1,CIiter+1)=CId(1,CIiter+1)-... CIa(1,CIiter)*CId(1,CIiter)/... CIb(1,CIiter); end for CIiter=1:CILen-2 if CIiter==1 CIf2(1, CILen-CIiter)=... CId(1,CILen-CIiter-1)/... CIb(1,CILen-CIiter-1); else CIf2(1, CILen-CIiter)=... (CId(1,CILen-CIiter-1)-... CIc(1,CILen-CIiter-1)*... CIf2(1, CILen-CIiter))/... CIb(1,CILen-CIiter-1); end end end %---------------- Generate CIf2 -------------------------- switch CItype case 1 CIf2(1,1)=3/CIh(1,1)*(CID(1,1)-... CIm1)-CIf2(1,2)/2; CIf2(1,CILen)=3/CIh(1,CILen-1)*... (CImN-CID(1,CILen-1))-... CIf2(1,CILen-1)/2; case 2 CIf2(1,1)=0.0; CIf2(1,CILen)=0.0; case 3 CIf2(1,1)=CIf2(1,1)-CIh(1,1)*(CIf2(1,3)-... CIf2(1,2))/CIh(1,2); CIf2(1,CILen)=CIf2(1,CILen-1)-... CIh(1,CILen-1)*(CIf2(1,CILen-1)-... CIf2(1,CILen-2))/CIh(1,CILen-2); case 4 CIf2(1,1)=CIf2(1,2); CIf2(1,CILen)=CIf2(1,CILen-1); case 5 CIf2(1,1)=CIm1; CIf2(1,CILen)=CImN; case 6 CIf2(1,1)=CIf2(1,CILen-1); CIf2(1,CILen)=CIf2(1,2); end if CIarr==0; CIcoes=zeros(CILen-1,4); for CIiter=1:CILen-1 CIcoes(CIiter,1)=CIY(CIiter,1); CIcoes(CIiter,2)=CID(1,CIiter)-CIh(1,CIiter)*... (2*CIf2(1,CIiter)+... CIf2(1,CIiter+1))/6; CIcoes(CIiter,3)=CIf2(1,CIiter)/2; CIcoes(CIiter,4)=(CIf2(1,CIiter+1)-... CIf2(1,CIiter))/... (6*CIh(1,CIiter)); end elseif CIarr==1 CIcoes=zeros(4,CILen-1); for CIiter=1:CILen-1 CIcoes(1,CIiter)=CIY(CIiter,1); CIcoes(2,CIiter)=CID(1,CIiter)-CIh(1,CIiter)*... (2*CIf2(1,CIiter)+CIf2(1,CIiter+1))/6; CIcoes(3,CIiter)=CIf2(1,CIiter)/2; CIcoes(4,CIiter)=(CIf2(1,CIiter+1)-... CIf2(1,CIiter))/(6*CIh(1,CIiter)); end end -- -- ~我不僅喜歡酒味嗆人的菸味更喜歡女人身體的耐人尋味 ~ -- ※ 發信站: 批踢踢實業坊(ptt.cc) ◆ From: 59.104.72.77 ※ 編輯: hilosi 來自: 59.104.72.77 (11/22 21:27) ※ 編輯: hilosi 來自: 59.104.72.77 (11/22 21:31)