How to solve a ode within a for loop?

49 vues (au cours des 30 derniers jours)
mathru
mathru le 10 Juil 2021
Commenté : mathru le 10 Juil 2021
I am trying to solve a second order diferential function using ode45. I have defined a function for the differential function where in the expresion there is a constant, say, A. For single value, code is working good. I have written my code as follows:
global A
A = 3;
tspan = 0: 0.1: 1;
r = 0.7;
r0 = 2;
c = [r r0];
[t1 p1] = ode45('height',tspan,c);
where the function is as follows:
global A
function pdot = height(t,p)
pdot = zeros(size(p))
pdot(1)=p(2);
B=5*A*p(1);
pdot(2)=B/5*p(2)^2-3*p(2)+6;
end
Now I want to solve the ode for different values of A using a for loop as follows:
global A
A = 0: 5: 10;
p32 = zeros(length(tspan),length(A));
for k = 1: length(A)
tspan = 0: 0.1: 1;
r = 0.7;
r0 = 2;
c = [r r0];
[t1 p1] = ode45('height',tspan,c);
p32(k) = p(:,1);
end
I am getting the following errors:
"Unable to perform assignment because the left and right sides have a different number of elements."
"f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0"
How can I solve the problem?

Réponse acceptée

Alan Stevens
Alan Stevens le 10 Juil 2021
Like this
A = 0:0.05:0.35;
tspan = 0: 0.1: 1;
p32 = zeros(length(tspan),length(A));
r = 0.7;
r0 = 2;
p0 = [r; r0];
for k = 1: length(A)
[t1, p1] = ode45(@(t,p) height(t,p,A(k)),tspan,p0);
p32(:,k) = p1(:,1);
end
plot(t1,p32),grid
xlabel('time'), ylabel('p32')
function pdot = height(~,p,A)
pdot = zeros(size(p));
pdot(1)=p(2);
B=5*A*p(1);
pdot(2)=B/5*p(2)^2-3*p(2)+6;
end
However, it fails with A larger than about 0.35. Have you got your equations right? (For example, should p(2)^2 be in the denominator?).
  1 commentaire
mathru
mathru le 10 Juil 2021
Hi Alan,
It works! Thanks for the code.
My original equation is too big. I just provided the first 3 terms.

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Loops and Conditional Statements dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by