plotting differential equation converging to stability value

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

Answer 1

[Star Strider]

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

[Joris]

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]

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.