How to solve a ode within a for loop?

54 Ansichten (letzte 30 Tage)
mathru
mathru am 10 Jul. 2021
Kommentiert: mathru am 10 Jul. 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?

Akzeptierte Antwort

Alan Stevens
Alan Stevens am 10 Jul. 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 Kommentar
mathru
mathru am 10 Jul. 2021
Hi Alan,
It works! Thanks for the code.
My original equation is too big. I just provided the first 3 terms.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Loops and Conditional Statements finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by