有人可以帮我解决这个代码吗?

问题描述 投票:0回答:1

我试图使用 interp1 引入一个变量 p 作为时间的函数,但它不起作用,它给了我一个向量,而不是得到一个值作为积分的结果,

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
matlab plot numerical-integration
1个回答
1
投票

您在使用

interp1
时遇到问题。

p = interp1(t,p,T)
- 这返回一个向量

'p = interp1(t,p,method,'pp')
- 这返回一个具有分段函数的结构。

文档在这里。

此外,您还可以使用

ppval
从分段函数中获取值。 最后一点 - Matlab r2014a 中有一个警告:

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

所以准备好和

GRIDDEDINTERPOLANT
一起去吧。

© www.soinside.com 2019 - 2024. All rights reserved.