現在この質問をフォロー中です
- フォローしているコンテンツ フィードに更新が表示されます。
- コミュニケーション基本設定に応じて電子メールを受け取ることができます。
How can I plot the derivatives of the components of the solution to a system of ODEs?
1 回表示 (過去 30 日間)
古いコメントを表示
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
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
2021 年 9 月 20 日
As always, my pleasure!
The code I posted runs without error. It should work with other differential equation functions as well, with changes in the loop appropriate to the function.
I only see the error message, not the code that produced it, so I have no idea what the problem is. If you copy and paste my code exactly as I wrote it, it should work correctly when you run it.
.
Cris19
2021 年 9 月 20 日
I just copied and pasted the commands you have written. On the interval [0,200] there is no error message, but I get that on the time interval [0,50]. I do not know how to fix it.
Star Strider
2021 年 9 月 20 日
Both of them work for me —
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', 'LineWidth',1.25)
plot(t, dwdt(3,:), ':b', 'LineWidth',1.25)
hold off
legend('$x(t)$','$y(t)$','$\dot{x}$','$\dot{y}$', 'Interpreter','latex', 'Location','best');
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/744469/image.png)
tspan = [0 200];
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', 'LineWidth',1.25)
plot(t, dwdt(3,:), ':b', 'LineWidth',1.25)
hold off
legend('$x(t)$','$y(t)$','$\dot{x}$','$\dot{y}$', 'Interpreter','latex', 'Location','best');
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/744474/image.png)
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
I made a few cosmetic changes to the plots to make the derivatives easier to see, and nothing else except to duplicate the code with the different values for ‘tspan’ in each example.
.
Cris19
2021 年 9 月 20 日
Thank you very much! It seems to work well for me on [0,50], but after... restarting the computer! Very strange...
Cris19
2021 年 9 月 23 日
編集済み: Cris19
2021 年 9 月 23 日
I would have one more question, if possible...
I tried to change f1 and f2 in the function 'coupled' with some piecewise functions:
function dwdt=complicate(t,w)
beta=1+exp(-t);
delta=2+exp(-t);
gama=exp(-t);
if t >= 1
f1 = 1/t;
else
f1 = -2*t^2+3*t;
end;
if t >= 2
f2 = 1/(t-1);
else
f2 = -(3/4)*t^2+2*t;
end;
h1=diff(f1(t),t);
h2=diff(f2(t),t);
r1=1/(t+2)^2;
r2=2/(t+1)^2;
dwdt=zeros(4,1);
dwdt(1)=w(2)-f1*w(1);
dwdt(2)=(h1+f1^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)+(h2+f2^2-delta)*w(3)-f2*w(4)-r2*w(3)^2;
end
where h1 and h2 are the derivatives of f1, f2. With the commands
tspan = [0 100];
z0=[0.01 0.01 0.01 0.01];
[t,z] = ode45(@(t,z) complicate(t,z), tspan, z0);
for k = 1:numel(t)
dwdt(:,k) = complicate(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');
I get the following errors
Array indices must be positive integers or logical values.
Error in complicate (line 15)
h1=diff(f1(t),t);
Error in @(t,z)complicate(t,z)
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
I do not know how can I fix these. It seems that not even z(:,1), z(:,3) can be plotted...
Star Strider
2021 年 9 月 23 日
With all due respect, there is so much that is wrong with the ‘complicate’ function that it will take a bit to describe them.
The direct answer is that it will just not work. At least not in its current form.
First, creating abrupt transitions creates non-differentiable ‘steps’ that no numerical differential equation integrator is able to deal with. The solution to that is to use an events function (see ODE Event Location) to stop the integration and re-start it with the previous step solution as the new initial conditions, and continue from there. Repeat for each additional ‘step’.
Second, the diff function (when used outside the Symbolic Math Toolbox, note the symbolic diff function) does not take the derivative, instead calculating the differences between adjacent elements of its argument vector. So here, it interprets ‘f(t)’ as its argument vector, ‘t’ as a subscript to it, and the second argument ‘t’ as the difference order.
So, ‘duplicate’ needs to become three different functions, one for t<1, one for t>=1 & t<2, and one for t>=2. Also, ‘h1’ and ‘h2’ need to be expressed appropriately, either by hand calculation or using the Symbolic Math Toolbox to differentiate them.
Correct these problems (I don’t see any others), use an event function, and it should work.
.
Cris19
2021 年 9 月 23 日
編集済み: Cris19
2021 年 9 月 23 日
Thanks a lot ! It is so complicated to me the using of event functions, even more so than the function 'complicate' itself... But I think it worked by calculating separately h1 and h2 as piecewise functions, replacing the lines
h1=diff(f1(t),t);
h2=diff(f2(t),t);
with
if t>=1
h1 = -1/t^2;
else
h1=-4*t+3;
end;
if t>=2
h2 = -1/(t-1)^2;
else
h2 = -(3/2)*t+2;
end;
I get the following plottings:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/747699/image.jpeg)
Can these plottings be considered correct?
Star Strider
2021 年 9 月 23 日
‘Can these plottings be considered correct?’
I’m not certain. The step changes in the function occur in the very beginning of the solution, so it’s difficult to see any effect. Most of this is after t=2, so there are no more ‘step’ discontinuities. The correct procedure is still to use events, however if there’s not much difference in the functions across the discontinuity, it may not be a problem.
.
Cris19
2021 年 9 月 23 日
Thank you very much for your answer. Do you refer to the discontinuities of h1, h2, because they are continuous on [0,infinity) ?
I would be very grateful if you could help me with the code using event functions, because, sincerely, I have never used these functions.
Star Strider
2021 年 9 月 23 日
‘Do you refer to the discontinuities of h1, h2, because they are continuous on [0,infinity) ?’
No. It is because they’re discontinuous at t=1 and t=2.
The event functions, although theoretically necessary if the discontinuities depend on the results of the integration, not time, may not be necessary here because I see no significant discontinuities in the results ar those points. In reality, since the values of the results do not depend on discontinuities in the results, only time, stopping the integration at t=1 and re-starting it using the results at t=1 and integrating to t=2 and doing the same there and then integrating through the rest of the vector may be all that is necessary, and no event functions would be needed. (Going back and studying the function and what it does, produce necessary insights such as these!)
.
Cris19
2021 年 9 月 23 日
編集済み: Cris19
2021 年 9 月 23 日
Ok, I try to understand... The left-hand limit of h1 at t=1 equals the right-hand limit at t=1 and equals -1; the same for h2 at t=2. So h1 and h2 are continuous at t=1, respectively t=2.
I would go further with h1 and h2 defined piecewisely. In my opinion, the plottings seem fine.
Star Strider
2021 年 9 月 23 日
I agree that it does not appear to be a problem. However if there are significant discontinuities in the functions across the discontinuities, stopping and re-starting would be necessary. I do not see that ‘f1’ and ‘f2’ or ‘h1’ and ‘h2’ are defined for t<1. I may be missing something since the code appears to work.
.
Cris19
2021 年 9 月 23 日
I have defined the new functions f1, f2
if t >= 1
f1 = 1/t;
else
f1 = -2*t^2+3*t;
end;
if t >= 2
f2 = 1/(t-1);
else
f2 = -(3/4)*t^2+2*t;
end;
and h1, h2
if t>=1
h1 = -1/t^2;
else
h1=-4*t+3;
end;
if t>=2
h2 = -1/(t-1)^2;
else
h2 = -(3/2)*t+2;
end;
in the new example, given in today's question.
I get no errors and the compiler does not seem to need restarting.
Thank you so very much, I have learned many new and interesting things and I hope I did not bother you too much.
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!エラーが発生しました
ページに変更が加えられたため、アクションを完了できません。ページを再度読み込み、更新された状態を確認してください。
Web サイトの選択
Web サイトを選択すると、翻訳されたコンテンツにアクセスし、地域のイベントやサービスを確認できます。現在の位置情報に基づき、次のサイトの選択を推奨します:
また、以下のリストから Web サイトを選択することもできます。
最適なサイトパフォーマンスの取得方法
中国のサイト (中国語または英語) を選択することで、最適なサイトパフォーマンスが得られます。その他の国の MathWorks のサイトは、お客様の地域からのアクセスが最適化されていません。
南北アメリカ
- América Latina (Español)
- Canada (English)
- United States (English)
ヨーロッパ
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
アジア太平洋地域
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)