MATLAB Answers

How can I plot the derivatives of the components of the solution to a system of ODEs?

조회 수: 60(최근 30일)
Cris19
Cris19 2021년 9월 20일
댓글: Star Strider 2021년 9월 23일 19:51
I have the following code for the function solving a system of ODEs:
function dwdt=coupled(t,w)
beta=1+exp(-t);
delta=2+exp(-t);
gama=exp(-t);
f1=1/(t+1);
f2=1/(t+2);
r1=1/(t+2)^2;
r2=2/(t+1)^2;
dwdt=zeros(4,1);
dwdt(1)=w(2)-f1*w(1);
dwdt(2)=-beta*w(1)-f1*w(2)+gama*w(3)-r1*w(1)^2;
dwdt(3)=w(4)-f2*w(3);
dwdt(4)=gama*w(1)-delta*w(3)-f2*w(4)-r2*w(3)^2;
xdot=w(2)-f1.*w(1);
end
With the commands
>> tspan = [0 50];
z0=[0.01 0.01 0.01 0.01];
[t,z] = ode45(@(t,z) coupled(t,z), tspan, z0);
figure
plot(t,z(:,1),'r',t,z(:,3),'b');
legend('x(t)','y(t)');
I have plotted the first w(1) and the third w(3) components of the solution. But I want also to plot on the same interval, with the same initial data, the derivatives of w(1) and w(3), i.e.,
dwdt(1)=w(2)-f1*w(1)
and
dwdt(3)=w(4)-f2*w(3)
If I add the command line
plot(t,z(:,1),'r',t,z(:,2)-f1*z(:,1),'b');
I get the message
Unrecognized function or variable 'f1'.
How dwdt(1) and dwdt(3) can be plotted ?

채택된 답변

Star Strider
Star Strider 2021년 9월 20일
The only way is to use al loop —
tspan = [0 50];
z0=[0.01 0.01 0.01 0.01];
[t,z] = ode45(@(t,z) coupled(t,z), tspan, z0);
for k = 1:numel(t)
dwdt(:,k) = coupled(t(k),z(k,:));
end
figure
plot(t,z(:,1),'r',t,z(:,3),'b')
hold on
plot(t, dwdt(1,:), ':r')
plot(t, dwdt(3,:), ':b')
hold off
legend('$x(t)$','$y(t)$','$\frac{dx(t)}{dt}$','$\frac{dy(t)}{dt}$', 'Interpreter','latex', 'Location','best');
function dwdt=coupled(t,w)
beta=1+exp(-t);
delta=2+exp(-t);
gama=exp(-t);
f1=1/(t+1);
f2=1/(t+2);
r1=1/(t+2)^2;
r2=2/(t+1)^2;
dwdt=zeros(4,1);
dwdt(1)=w(2)-f1*w(1);
dwdt(2)=-beta*w(1)-f1*w(2)+gama*w(3)-r1*w(1)^2;
dwdt(3)=w(4)-f2*w(3);
dwdt(4)=gama*w(1)-delta*w(3)-f2*w(4)-r2*w(3)^2;
xdot=w(2)-f1.*w(1);
end
Experiment to get different results.
.
  댓글 수: 20
Star Strider
Star Strider 2021년 9월 23일 19:51
As always, my pleasure!
No worries, no bother!
.

댓글을 달려면 로그인하십시오.

추가 답변(0개)

태그

Community Treasure Hunt

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

Start Hunting!

Translated by