フィルターのクリア

How to solve a ode within a for loop?

18 ビュー (過去 30 日間)
mathru
mathru 2021 年 7 月 10 日
コメント済み: mathru 2021 年 7 月 10 日
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?

採用された回答

Alan Stevens
Alan Stevens 2021 年 7 月 10 日
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 件のコメント
mathru
mathru 2021 年 7 月 10 日
Hi Alan,
It works! Thanks for the code.
My original equation is too big. I just provided the first 3 terms.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeLoops and Conditional Statements についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by