Может ли кто-нибудь помочь мне решить этот код, пожалуйста?

Question

Я пытаюсь ввести переменную р как функцию времени с помощью interp1, но он не работает, и это дает мне вектор вместо того, чтобы одно значение как результат для моей интеграции,

function RunLogisticOscilFisher 

omega = 1; 

k = 10; 

N0 = 1; 

A = 1; 

p0 = .1; 

tspan=(0:0.1:10); 

% Finding the numerical solution for the function using ode45 solver 

[t,p]=ode45(@logisticOscilfisher,tspan,p0,[],N0,k,omega); 

% [t,p]=ode23s(@(t,p) N0*sin(omega*t)*p*(1-p./k),tspan,p0,odeset('AbsTol',1e-8,'RelTol',1e-10')); 

% Plotting the function with time 

figure(1) 

plot(t,p) 

% Finding the integral to get the fisher information with respect to p 

f = @(p) ((A.*(((N0*sin(omega*t).^2.*(1-2*p./k))+(omega.*cos(omega*t)) 

    ).^2)./(N0.^2*sin(omega*t).^4.*(p-p.^2./k).^2))) 

I1=integral(f, 11,20,'ArrayValued',true) 

I2=integral(f,11,40,'ArrayValued',true) 

I3=integral(f,11,60,'ArrayValued',true) 

I4=integral(f,11,80,'ArrayValued',true) 

I5=integral(f,11,101,'ArrayValued',true) 

I=[I1,I2,I3,I4,I5] 

II=[I1./20 I2./40 I3./60 I4./80 I5./101] 

T=[20 40 60 80 101]'; 

%Plotting the Fisher Information 

figure(2) 

plot(T,I,'b-'), hold on 

plot(T,II,'r-') 

hold off 

% Finding the integral to get the fisher information with respect to t 

P = @(T) interp1(t,p,T); 

ff = @(t) (A.*(((N0*sin(omega*t).^2.*(1-2*p./k))+(omega.*cos(omega*t)) 

     ).^2)./(N0.^2*sin(omega*t).^4.*(p-p.^2./k).^2)) 

J1=integral(ff, 0.1,2,'ArrayValued',true) 

J2=integral(ff, 0.1,4,'ArrayValued',true) 

J3=integral(ff, 0.1,6,'ArrayValued',true) 

J4=integral(ff,0.1,8,'ArrayValued',true) 

J5=integral(ff,0.1,10,'ArrayValued',true) 

J=[J1,J2,J3,J4,J5] 

JJ=[J1./2 J2./4 J3./6 J4./8 J5./10] 

R=[2,4,6,8,10]'; 

%Plotting the Fisher Information 

figure(3) 

plot(R,J,'b-'), hold on 

plot(R,JJ,'r-') 

hold off 

figure(4) 

plot(t,f(t)) 

figure(5) 

plot(log(t),log(f(t))) 

P = @(T) interp1(t,p,T); 

a=exp(1-cos(t)); 

fff = @(T) ( (A.* ((sin(t)).^2 .* (1-2.*(a./ (9.9 + 0.1 .* a))./10) + cos(t)).^2 

     ) ./ ((sin(t)).^4 .* ((a ./ (9.9+(0.1.*a))) - ((a ./ (9.9 + (0.1 .* a)) 

     ).^2 ./10)).^2 ) ) 


K1=integral(fff, 0,2,'ArrayValued',true) 

K2=integral(fff, 0,4,'ArrayValued',true) 

K3=integral(fff, 0,6,'ArrayValued',true) 

K4=integral(fff,0,8,'ArrayValued',true) 

K5=integral(fff,0,10,'ArrayValued',true) 

K=[K1,K2,K3,K4,K5] 

KK=[K1./2 K2./4 K3./6 K4./8 K5./10] 

%Plotting the Fisher Information 

figure(6) 

plot(R,K,'b-'), hold on 

plot(R,KK,'r-') 

hold off 

1; 

% function dpdt = logisticOscilfisher(t,p,N0,k,omega) 

% dpdt = N0*sin(omega*t)*p*(1-p/k); 

% end 

Answer 1

У вас есть проблема с использованием interp1.

p = interp1(t,p,T) - это возвращает вектор 'p = interp1(t,p,method,'pp') - это возвращает структуру с кусочно-мутной функцией.

Docs находятся здесь.

Кроме того, вы можете использовать ppval для получения значений от кусочной функции. И последнее замечание - есть предупреждение в Matlab r2014a:

Warning: INTERP1(...,'PP') will be removed in a future release. Use GRIDDEDINTERPOLANT instead. 

Так что будьте готовы идти с GRIDDEDINTERPOLANT.