plotting differential equation converging to stability value
5 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/168305/image.jpeg)
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
0 commentaires
Réponses (1)
Star Strider
le 30 Déc 2017
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 commentaires
Star Strider
le 30 Déc 2017
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.
Voir également
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!