How to solve a ode within a for loop?

47 views (last 30 days)
mathru
mathru on 10 Jul 2021
Commented: mathru on 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?

Accepted Answer

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

Sign in to comment.

More Answers (0)

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!