Как решить ODE серии, которые включают суммирование - PullRequest
Новая
       3

Как решить ODE серии, которые включают суммирование

0 голосов
/ 05 марта 2020

Я пытаюсь решить систему ODE в Matlab. Я пробовал несколько способов включить сумму суммирования, но понятия не имею, как это сделать. Я получаю следующую ошибку.

Неопределенная функция 'symsum' для входных аргументов типа 'double'.

мой код: Прикреплено фактическое уравнение модели 

imax=0;
syms k
m=10; %mosquito per human
a=.21; %bitting frequiency (per day)
t0=0;
tmax=100;
g=1/10;      %mosquito death rate (per day)
n=12;        %duration of sporogony in mopsquito
b=.5;        %transmission probability: mosquito to human 
c=.23;       %transmission probability: human to mosquito
alpha=1/332;  %hypnozoites activation rate
r=0;        %rate of blood stage infection clearance 
mu=1/425  ;   %hypnozoites death rate
Nval=8.5;
N=min(imax,Nval);         %Number of hypnozoites per infection

%initial condition 
S0=.7;
I0=.2;
S1=0;
I1=0;
S2=0;
I2=0;
Sm=.99;
Im=.01;


for i=1:imax+1

    options = odeset('RelTol', 1e-5);
% y0=[1-1e-6,1e-6,0];
    [t, y]=ode45(@vivax_model2,[t0,tmax],[S0,I0,Sm,Im],options,[m a g n b c alpha mu N r]);
    plot(t,y(:,i),t,y(:,i+imax+1),t,y(:,3*imax+1),t,y(:,3*imax+2))
    hold on
    xlabel('Time (Days)')
    ylabel('Fraction of population in each Class')
%     legend('Susceptible human with no hypno','Susceptible human with 1 hypno','Susceptible human with 2 hypno','Infected human with no hypno','Infected human with 1 hypno','Infected human with 2 hypno','susceptible mosquito','Infected mosquito')

end

function ydot=vivax_model2(t,y,parameter)

m=parameter(1); a=parameter(2); g=parameter(3);n=parameter(4);b=parameter(5);
c=parameter(6);alpha=parameter(7);mu=parameter(8);N=parameter(9);r=parameter(10);
imax=0;
sum=0;
for i=1:imax+1
    ydot=zeros(2*(imax+1)+2,1);
    lambda=m*a*b*y(2*(imax+1)+2);
    ydot(i)=-lambda*y(i)-i*(mu+alpha)*y(i)+(i+1)*mu*y(i+1)+r*r*y(i+imax+1);
    for k=1:imax+1

        ydot(i+imax+1)=-lambda*y(i+imax+1)+symsum(lambda*((N/(N+1)).^(i-k))*(1/(N+1))*(y(k)+y(imax+1+k)),k,1,imax+1)...
                   -i*(mu+alpha)*y(i+imax+1)+(i+1)*(mu+alpha)*y(i+1+(imax+1))+(i+1)*alpha*y(i+1)-r*y(i+imax+1);
    end             
    ydot(2*(imax+1)+1)=g-a*c*symsum((lambda*(y(k+imax+1))),k,1,imax+1)*(exp(-g*n)-y(3*imax+2))-g*y(2*(imax+1)+1);
    ydot(2*(imax+1)+2)=a*c*symsum((lambda*(y(k+imax+1))),k,1,imax+1)*(exp(-g*n)-y(3*imax+2))-g*y(2*(imax+1)+2);

end
...