what's the problem in defining the function in for loop?
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
SAHIL SAHOO
am 31 Jul. 2022
Kommentiert: Walter Roberson
am 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 Kommentare
Akzeptierte Antwort
Walter Roberson
am 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 Kommentare
Walter Roberson
am 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')
Weitere Antworten (0)
Siehe auch
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!