How to remove points in time vector when I am getting a blowup Fisher information at those points please?

2 次查看（过去 30 天）

Avan Al-Saffar 2014-12-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/167095-how-to-remove-points-in-time-vector-when-i-am-getting-a-blowup-fisher-information-at-those-points-p

编辑： dpb 2014-12-18

The code is :

function RunlogOscilnumericfishfixedn0omega2
omega=1;
N0=1;
k = 10;
A = 1;
p0 = 0.1;
t=(0.01:0.1:1000);
n=length(t);
[t,p] = ode45(@logOscilnumeric2,t,p0,[],omega,k,N0);
P = @(T) interp1(t,p,T);
v = (N0.^2.*P(t).^2.*(k.*sin(omega.*t) - P(t)).^4);
for i = 1 : n
 if abs(v(i)) < 1e-15
    f_v(i) = 0;
 else
 f_v(i) = ( (A.*( (k.^2 .*N0.*sin(omega.*t(i))) - (3.*N0.*P(t(i)).*k.*sin(omega.*t(i))) + 
(k.^2.*omega.*cos(omega.*t(i))) + (2.*N0.*P(t(i)).^2) ).^2 )./v(i) ) ;
 end
end
f = @(t) f_v;
I = integral( f, 0.01,1000,'ArrayValued',true)./1000
figure(1)
plot(t,I)
xlabel('Time')
ylabel('Fisher Information')
P = @(T) interp1(t,p,T);
g = @(t) ( (A.*( (k.^2 .*N0.*sin(omega.*t)) - (3.*N0.*P(t).*k.*sin(omega.*t)) + 
(k.^2.*omega.*cos(omega.*t)) + (2.*N0.*P(t).^2) ).^2 )./ (N0.^2.*P(t).^2.*(k.*sin(omega.*t) - P
(t)).^4)   );
figure(2)
plot(t,g(t))
1;
% function dpdt = logOscilnumeric2(t,p,omega,k,N0)
% dpdt = N0*sin(omega*t)*p - (N0 *p.^2/k);
% end