plotting differential equation converging to stability value

3 次查看（过去 30 天）
Joris
Joris 2017-12-30
评论： Star Strider 2017-12-30
Hello everyone, I have a question about plotting a differential equation converging to a stability value. I have the following equation:
dpdt = (alpha*cos(phi(ii))-beta*sin(phi(ii))-1)/(1-cos(phi(ii)))
Given a certain value for alpha and beta, the value of phi converges to a steady state equilibrium. To determine the way it goes to this value I've created the following script:
clear all; close all; clc;
a = [0.3 0.35 0.4 0.45 0.5 0.6];
b = 1.1;
for i = 1:6
alpha = a(i);
x = [1:10000001];
dpdt = zeros(1,10000000);
phi = zeros(1,1000001);
phi(1) = pi/4;
for ii = 1:10000000
dpdt(ii) = (alpha*cos(phi(ii))-b*sin(phi(ii))-1)/(1-cos(phi(ii)));
phi(ii+1) = phi(ii) + dpdt(ii)/10000000;
end
plot(x/10000000,phi)
xlabel('angle theta [rad]')
ylabel('angle phi [rad]')
hold on
end
Now the problem with this script is that it has a certain accuracy based on the value of x (1000000 or higher etc.) This seems to me that it is not a very efficient way of plotting the differential equation. I've looked into using the function 'dsolve' but I was not able to get this to work. The plot that I'm getting now (takes about 1 minute) looks like this:
In which you can see that the accuracy is rather low since the value of the slope for the lower lines. An extra zero seems like an inefficient solution. This has got to be easier right? Does anyone have any solutions?
thanks in advance, Joris

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

回答（1 个）

Star Strider
Star Strider 2017-12-30
Use the efficient and robust core MATLAB differential equation solvers instead of your loop.
I am not certain that I coded your differential equation correctly, so make any necessary changes to it.
See if this does what you want:
dpdt = @(t,phi,alpha,beta) [phi(1); (alpha.*cos(phi(1))-beta*sin(phi(1))-1)/(1-cos(phi(1)))];
a = [0.3 0.35 0.4 0.45 0.5 0.6];
b = 1.1;
tspan = linspace(0, 2, 50);
for k1 = 1:length(a)
[t,thphi{k1}] = ode45(@(t,phi) dpdt(t,phi,a(k1),b), tspan, [pi/4; -0.25]);
end
AxH = axes('nextplot', 'add');
figure(1)
for k1 = 1:numel(a);
plot(AxH, thphi{k1}(:,1), thphi{k1}(:,2))
end
grid
  2 个评论
Joris
Joris 2017-12-30
编辑：Joris 2017-12-30
Thanks! this looks very promising. There are some problems with it however. To explain let me tell you more about the problem at hand. The 't' is an angle 'theta' so the differential equation describes 'phi' differentiated w.r.t. 'theta' The problem looks like this:
Where the red lines represent a socalled log-truck and A is the front of the truck (carrying logs through a corner). Hopefully its more clear now. After adjusting the script a bit to fit the initial conditions better, this is the result:
Now, the problem is that this convergence (if there is any) looks really weird. One would expect that when the truck is driving along on the roundabout a stable value of phi will occur. This does not occur right now. I'm currently looking further into it but now, with the additional information, do you have any suggestions? Although your script works a lot more fluid, I would expect the results to be more of my plot in the initial question.
This is how I implemented the script right now:
dpdt = @(theta,phi,alpha,beta) [phi(1); (alpha.*cos(phi(1))-beta*sin(phi(1))-1)/(1-cos(phi(1)))];
a = [0.3 0.35 0.4 0.45 0.5 0.6];
b = 2;
thetaspan = linspace(0, 2, 50);
for k1 = 1:length(a)
[theta,thphi{k1}] = ode45(@(theta,phi) dpdt(theta,phi,a(k1),b), thetaspan, [pi/2; pi/4]);
end
AxH = axes('nextplot', 'add');
figure(1)
for k1 = 1:numel(a)
plot(AxH, thphi{k1}(:,1), thphi{k1}(:,2))
end
grid
Thanks a lot for the quick reply though!
Star Strider
Star Strider 2017-12-30
My pleasure.
I have no idea what you are doing. My background is biomedical engineering, and while I would like to help, the mechanical engineering aspects of your problem are significantly outside my areas of expertise.
I will of course help you with the MATLAB code and integrating your differential equations, however I cannot help with respect to the engineering, and coding your differential equations.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by