what's the problem in defining the function in for loop?
1 view (last 30 days)
Show older comments
SAHIL SAHOO
on 31 Jul 2022
Commented: Walter Roberson
on 31 Jul 2022
clc
ti = 0; %inital time
tf = 70E-9;% final time
tspan=[ti tf];
tp =1E-12;
T = 2E3;
p = 0.05;
k = (0.62).*10^(4);
c = 3E8;
N = 3.0;
n = k*(c/N)*tp
a = 5;
f = @(t,y) [
(y(4).*y(1) - n.*y(2).*sin(y(3)));
(y(5).*y(2) + n.*y(1).*sin(y(3)));
(-a.*(y(5)-y(4)) + n.*((y(1)./y(2)) - (y(2)./y(1))).*cos(y(3)));
(p - y(4)-(1+2.*y(4)).*(y(1)^2))./T;
(p - y(5)-(1+2.*y(5)).*(y(2)^2))./T;
];
[time,Y] = ode45(f,tspan./tp,[sqrt(p);sqrt(p);0;0;0.01]);
O = linspace(-10,10,100);
U = zeros(length(O),1) ;
for i = 1:length(O)
U(i) = -a.*(Y(:,5) - Y(:,4)).*O(i) + n.*((Y(:,1)./Y(:,2)) - (Y(:,2)./Y(:,1)))*sin(O(i));
end
subplot 311
plot(time,Y(:,3));
xlabel('time')
ylabel('phase')
subplot 312
plot(time,Y(:,2));
xlabel('time')
ylabel('Amplitude')
subplot 312
plot(O,U);
xlabel('phase')
ylabel('potential')
0 Comments
Accepted Answer
Walter Roberson
on 31 Jul 2022
tspan=[ti tf];
[time,Y] = ode45(f,tspan./tp,[sqrt(p);sqrt(p);0;0;0.01]);
When you use a vector of length 2 for the tspan, ode45 generates time outputs whenever it feels like it, so time will be a vector whose length is not easily predictable ahead of time. The Y output will have as many rows as there are entries in time and will have as many columns as there are entries in the initial conditions.
U(i) = -a.*(Y(:,5) - Y(:,4)).*O(i) + n.*((Y(:,1)./Y(:,2)) - (Y(:,2)./Y(:,1)))*sin(O(i));
The left hand side, U(i), is a scalar location.
On the right hand side, you use Y(:,1), Y(:,2), Y(:,4), and Y(:,5), each of which is a column vector with as many entries as the number of time values that were returned by ode45(). So the right hand side is a column vector of results, and you are trying to fit the column vector into a scalar location.
3 Comments
Walter Roberson
on 31 Jul 2022
ti = 0; %inital time
tf = 70E-9;% final time
tspan = linspace(ti, tf, 1000);
%tspan = [ti, tf];
tp =1E-12;
T = 2E3;
p = 0.05;
k = (0.62).*10^(4);
c = 3E8;
N = 3.0;
n = k*(c/N)*tp
a = 5;
f = @(t,y) [
(y(4).*y(1) - n.*y(2).*sin(y(3)));
(y(5).*y(2) + n.*y(1).*sin(y(3)));
(-a.*(y(5)-y(4)) + n.*((y(1)./y(2)) - (y(2)./y(1))).*cos(y(3)));
(p - y(4)-(1+2.*y(4)).*(y(1)^2))./T;
(p - y(5)-(1+2.*y(5)).*(y(2)^2))./T;
];
[time,Y] = ode45(f,tspan./tp,[sqrt(p);sqrt(p);0;0;0.01]);
O = linspace(-0.08,0.08,10);
U = -a.*(Y(:,5) - Y(:,4)).*O + n.*((Y(:,1)./Y(:,2)) - (Y(:,2)./Y(:,1))).*sin(O);
subplot 311
plot(time,Y(:,3));
xlabel('time')
ylabel('phase')
subplot 312
plot(time,Y(:,2));
xlabel('time')
ylabel('Amplitude')
subplot 313
plot(O,U);
xlabel('phase')
ylabel('potential')
More Answers (0)
See Also
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!