2015-05-10 2 views
0

Мой сценарий должен запускать Runge-Kutta, а затем интерполировать вокруг вершин с помощью polyfit для вычисления максимальных значений вершин. Похоже, что значения х-значений максимальных точек правильны, но значения у-значения по какой-то причине отключены. Сидели с ним уже 3 дня. Проблема должна быть в последнем for-loop, когда я вычисляю py?Интерполяция с использованием 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) 

ответ

0

Ваша polyfit строка должна гласить:

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

Более интересно, принять к сведению предупреждения вы получите о бедные (я полагаю, вы видите это). Зачем? Частично из-за числовой точности (ваши цифры очень маленькие, масштабируются около 10^-6), а отчасти потому, что вы просите 4-й порядок соответствовать пяти пунктам (что единственно). Чтобы сделать это «лучше», используйте больше входных точек (более 5) или полиномиальную посадку более низкого порядка (квадратичная, как правило, много) и (возможно) перемасштабируйте, прежде чем использовать инструмент «полифит».

Сказав, что на практике эта проблема часто решается с использованием трех точек и квадратичной подгонки, поскольку она является вычислительно дешевой и дает почти те же ответы, что и более сложные подходы, но вы не получили этого от меня (с такие бесшумные данные, это не имеет большого значения в любом случае).

Смежные вопросы