使用polyfit插值（Matlab）

我的腳本應該運行Runge-Kutta，然後使用polyfit插入頂部以計算頂部的最大值。我似乎得到了最大點的x值正確，但y值由於某種原因而關閉。已經坐了3天了。問題應該是在我計算py的最後一個for循環中？使用polyfit插值（Matlab）

功能：

function funk = FU(t,u) 

L0 = 1; 
C = 1*10^-6; 

funk = [u(2); 2.*u(1).*u(2).^2./(1+u(1).^2) - u(1).*(1+u(1).^2)./(L0.*C)];

計劃：

%Runge kutta 

clear all 
close all 
clc 
clf 

%Given values 
U0 = [240 1200 2400]; 
L0 = 1; 
C = 1*10^-6; 

T = 0.003; 
h = 0.000001; 

W = []; 

% Runge-Kutta 4 

for i = 1:3 

    u0 = [0;U0(i)]; 

    u = u0; 
    U = u; 

    tt = 0:h:T; 

    for t=tt(1:end-1) 

     k1 = FU(t,u); 
     k2 = FU(t+0.5*h,u+0.5*h*k1); 
     k3 = FU((t+0.5*h),(u+0.5*h*k2)); 
     k4 = FU((t+h),(u+k3*h)); 

     u = u + (1/6)*(k1+2*k2+2*k3+k4)*h; 

     U = [U u]; 
    end 

    W = [W;U]; 

end 

I1 = W(1,:); I2 = W(3,:); I3 = W(5,:); 
dI1 = W(2,:); dI2 = W(4,:); dI3 = W(6,:); 

I = [I1; I2; I3]; 
dI = [dI1; dI2; dI3]; 

%Plot of the currents 
figure (1) 
plot(tt,I1,'r',tt,I2,'b',tt,I3,'g') 
hold on 
legend('U0 = 240','U0 = 1200','U0 = 2400') 

BB = []; 
d = 2; 
px = []; 
py = []; 

format short 

for l = 1:3 

    [H,Index(l)]=max(I(l,:)); 

    Area=[(Index(l)-2:Index(l)+2)*h]; 
    p = polyfit(Area,I(Index(l)-2:Index(l)+2),4); 

    rotp(1,:) = roots([4*p(1),3*p(2),2*p(3),p(4)]); 

    B = rotp(1,2); 

    BB = [BB B]; 

    Imax(l,:)=p(1).*B.^4+p(2).*B.^3+p(3).*B.^2+p(4).*B+p(5); 

    Tsv(i)=4*rotp(1,l); 


    %px1 = linspace(h*(Index(l)-d-1),h*(Index(l)+d-2)); 
    px1 = BB; 
    py1 = polyval(p,px1(1,l)); 

    px = [px px1]; 
    py = [py py1]; 

end 

% Plots the max points 
    figure(1) 
    plot(px1(1),py(1),'b*-',px1(2),py(2),'b*-',px1(3),py(3),'b*-') 
    hold on 

disp(Imax)

來源

2015-05-10 stx1020

你polyfit行應爲：

p = polyfit(Area,I(l, Index(l)-2:Index(l)+2),4);

更有趣的是，帶你獲取有關警告的音符較差的該多項式的條件（我假設你看到這些）。爲什麼？部分原因是數值精度（你的數字非常小，約爲10^-6），部分原因是你要求四階擬合爲五點（這是單數）。要做到這一點「更好」，使用更多的輸入點（超過5），或更低階的多項式擬合（二次方通常很多），並且可能在使用polyfit工具之前重新調整。

話雖如此，實際上這個問題通常使用三點和二次擬合來解決，因爲它在計算上便宜並且提供了與更復雜方法幾乎相同的答案，但是你沒有從我那裏得到（與像這樣的無噪聲數據，反正並不重要）。

來源

2015-05-11 17:12:31

使用polyfit插值（Matlab）

回答

相關問題