Help me with code to solve for loop to solve ODE

Hello everybody,

I have a problem with my code. I want to write a code to solve ODE (image). I am trying it with function ode45, I can only write code when value s=const. But I do not how to write a code for s(t) in this case. Can you show me how to do in this case. Thanks so much! Best regards!

Answer 1

[Ameer Hamza]

Try this code

tspan = [0 500];
ic = [1; 0; 0; 0; 0];
[t, p] = ode45(@odeFun, tspan, ic);

plot(t, p)
legend({'p', 'p\_dot1', 'p\_dot2', 'p\_dot3', 'p\_dot4'})

function dpdt = odeFun(t, p)

A = [-1 0 0 0 0;
     1 -1 0 0 0;
     0 1 -1 0 0;
     0 0 1 -1 0;
     0 0 0 1 -1];

B = [0 1 0 0 0;
     0 -1 1 0 0;
     0 0 -1 1 0;
     0 0 0 -1 1;
     0 0 0 0 -1];

lambda = 0.2;
T = 100;
g = T/2;
s0 = 0.5;

if (t-floor(t/T)*T)<g
    s = s0;
else
    s = 0;
end

dpdt = (lambda*A + s*B)*p;

end

[Le Duc Long]

Thanks so much Ameer Hamza,

In your code, What would the graph of s (t) look like?

This is my code for s(t):

syms t

s = piecewise((t-floor(t/T)*T)<g,s0,0);

fplot(s, [0 500])

Please see my code for s(t). And Can you show me in general case, the number of differential equations is n?

Best regards!

[Ameer Hamza]

Yes, s(t) is the same in my code, as shown in the figure. You can use if-else statements to construct an arbitrary piecewise signal.

[Le Duc Long]

Thanks again Bro!

[Ameer Hamza]

I am glad to be of help!

[Le Duc Long]

Hi again Ameer Hamza,

I am new member in Mathwork so i do not have some skill. I can ask you one more question for this problem?

1. In your code: "legend({'p', 'p\_dot1', 'p\_dot2', 'p\_dot3', 'p\_dot4'})". This is value p1(t), p2(t)....? So I think that code "legend({'p\_1', 'p\_2', 'p\_3', 'p\_4', 'p\_5'})" will be clearer. Do you think so?
2. I tried to draw graphic S(t). But i do not how to write code? Can you show me?

[Ameer Hamza]

Yes, in that case, this notation seems better. You can write it like this too

legend({'p_1', 'p_2', 'p_3', 'p_4', 'p_5'})

So that the number appears in the subscript.

[Le Duc Long]

Yeah. Can you show me method for the second question?

[Le Duc Long]

I tried with the FOR loop. But i do not know how to discrible pk(t)?

[Le Duc Long]

Hi Ameer Hamza,

I wrote code for general case. Can you see it:

clear all
tspan = [0:1:450];
n=12;
ic=zeros(1,n+2);
for ii=1:n+2
    ic(1,n+2)=1;
end
[t, p] = ode45(@odeFun, tspan, ic);
figure ('name','xac suat theo time')
plot(t, p(:,1),t,p(:,length(ic)))
legend({'p0', 'p12'})

r=zeros(length(tspan),1);
for t=1:length(tspan)
    r(t,1)=12;
    for k=1:(length(ic)-2)
        r(t,1)=r(t,1)+k*p(t,k+2);
    end
end
figure ('name','chieu dai hang cho theo time')
t=tspan ;
plot(t,r);
legend ({'hang cho'})

function dpdt = odeFun(t, p)
n=12;
for ii=1:n+2
    for jj=1:n+1
        if ii==jj A1(ii,jj)=-1;
        elseif ii==jj+1 A1(ii,jj)=1;
        else A1(ii,jj)=0;
        end
    end
end
A0=zeros(n+2,1);
A=[A1,A0];
for ii=1:n+2
    for jj=1:n+1
        if ii==jj B1(ii,jj)=1;
        elseif ii==jj+1 B1(ii,jj)=-1;
        else B1(ii,jj)=0;
        end
    end
end
B0=zeros(n+2,1);
B=[B0,B1];

lambda = 0.2;
T = 50;
g = T/2;
s0 = 0.5;

if (t-floor(t/T)*T)<g
    s = 0;
else
    s = s0;
end

dpdt = (lambda.*A + s.*B)*p;

end

[Le Duc Long]

Ameer Hamza, you can tell me when value lamda is not const, lamda has distribute Poisson. What kind of fuction I can use? Thanks so much!

[Ameer Hamza]

Try using https://www.mathworks.com/help/stats/poissrnd.html to generate random number using poisson distribution.

[Le Duc Long]

Thanks you very much Ameer Hamza.