Help me with code to solve for loop to solve ODE

2020 6 23
1 个回答
回答已采纳
更新时间：2020 6 25
3 次查看（30 天）

0 个投票

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!

请先登录，再回答此问题。

 采纳的回答

Ameer Hamza
Ameer Hamza 2020-6-23
编辑：Ameer Hamza 2020-6-23

0 个投票

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

12 个评论
显示 10更早的评论 隐藏 10更早的评论

Le Duc Long
Le Duc Long 2020-6-23
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
Ameer Hamza 2020-6-24
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
Le Duc Long 2020-6-24
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
Ameer Hamza 2020-6-24
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
Le Duc Long 2020-6-24
Yeah. Can you show me method for the second question?
Le Duc Long
Le Duc Long 2020-6-24
I tried with the FOR loop. But i do not know how to discrible pk(t)?
Le Duc Long
Le Duc Long 2020-6-25
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
Le Duc Long 2020-6-25
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
Ameer Hamza 2020-6-25
Try using https://www.mathworks.com/help/stats/poissrnd.html to generate random number using poisson distribution.
Le Duc Long
Le Duc Long 2020-6-25
Thanks you very much Ameer Hamza.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by