Using ode45 and plot
显示 更早的评论
I am working on 1. I am having trouble coding the function so that I can use the ode45 function. This is what my code looks like :
function xp=F(t,x)
xp=zeros(2,1); % since output must be a column vector
xp(2)=-(x(1)+(0*(x(1)^3)));
xp(3)=-(x(1)+(0.1*(x(1)^3)));
xp(4)=-(x(1)+(0.2*(x(1)^3)));
xp(5)=-(x(1)+(0.4*(x(1)^3)));
xp(6)=-(x(1)+(0.6*(x(1)^3)));
xp(7)=-(x(1)+(0.8*(x(1)^3)));
xp(8)=-(x(1)+(1.0*(x(1)^3)));
then in the command window I am inputting :
[t,x]=ode45('F',[0,20],[0,1]); plot(t,x(:,1))
I know I probably messed up the code . Any ideas?
采纳的回答
Star Strider
2015-7-10
Don’t do everything at once!
Otherwise, you’re almost there! This is how I would do it:
F = @(t,x,e) [x(2); -(x(1)+(e.*(x(1).^3)))]; % Spring ODE
ev = 0:0.2:1.0; % ‘epsilon’ Vector
ts = linspace(0,10, 50); % Time Vector
for k1 = 1:length(ev)
[~, x{k1}]=ode45(@(t,x) F(t,x,ev(k1)), ts, [0,1]);
end
figure(1)
for k1 = 1:length(ev)
subplot(length(ev), 1, k1)
plot(ts,x{k1})
legend(sprintf('\\epsilon = %.1f', ev(k1)))
grid
axis([xlim -2 +2])
end
xlabel('Time')
I stored ‘x’ as a cell array because it’s easier than storing it as a multi-dimensional array. The first variable, ‘x(:,1)’ is blue and ‘x(:,2)’ is red in each plot. I used subplots because it’s easier to compare the plots that way. I also specified the time argument, ‘ts’ as a vector so that all the integrated values would be the same size. These choices are for convenience, and not required.
9 个评论
Sam
2015-7-10
Star Strider
2015-7-10
My pleasure!
I got similar results when I (just now) varied epsilon from 0 to 50 in steps of 10. The value of epsilon affects the frequency, from 1/6.25 for epsilon = 0 to about 1/2.25 for epsilon = 50, while the amplitudes of ‘x(:,1)’ decreased from ±1 to ±0.42 as epsilon increased, while ‘x(:,2)’ remained stable at about ±1.
So you are correct — the amplitude of ‘x(:,1)’ varies inversely with epsilon, and the frequency varies proportionally with epsilon.
Sam
2015-7-12
Star Strider
2015-7-12
编辑：Star Strider
2015-7-12
Here is the way I usually use to detect zero-crossings:
t = 10:12;
x = [2 0 -2]';
zx = t(x .* circshift(x,[0 1]) <= 0)
The ‘zx’ assignment here gives the ‘t’ value for the zero-crossing. The ‘x’ vector does not have to include zero (it did here for illustration purposes) but the ‘x’ variable must change signs.
This will give you the approximate location of the zero-crossing. To get a more precise value, use an interpolation routine. For a detailed discussion of the techinque, see Curve Fitting to a Sinusoidal Function.
Note: Column vectors and row vectors have to be matched to the circshift second argument.
———————————————————————————————————
EDIT — This version of my original code will find and plot the initial zero-crossings for x(:,1). (In the ‘ixrng’ assignment, using the third value of ‘zxix’ was necessary for the code to avoid the initial value as the first zero.)
F = @(t,x,e) [x(2); -(x(1)+(e.*(x(1).^3)))]; % Spring ODE
ev = 0:0.2:1.0; % ‘epsilon’ Vector
ev = [0:1:5]*10;
ts = linspace(0,6.5, 50); % Time Vector
for k1 = 1:length(ev)
[~, x{k1}]=ode45(@(t,x) F(t,x,ev(k1)), ts, [0,1]);
end
figure(1)
for k1 = 1:length(ev)
subplot(length(ev), 1, k1)
plot(ts,x{k1})
hold on
legend(sprintf('\\epsilon = %.1f', ev(k1)))
grid
axis([xlim -2 +2])
xv = x{k1}(:,1);
zxix = find((xv.*circshift(xv, [1 0])) <= 0); % Approximate Zero-Crossing Indices
ixrng = max(1,zxix(3)-1):min(zxix(3)+1,length(ts)); % Index Range
tzx{k1} = interp1(xv(ixrng), ts(ixrng), 0); % ‘ts’ Values At Zero-Crossings
plot(tzx{k1}, zeros(size(tzx{k1})), 'bp') % Plot Zero Crossings
hold off
end
xlabel('Time')
Sam
2015-7-12
Star Strider
2015-7-12
My pleasure!
For the particular solution, the function changes to:
Fp = @(t,x,w) [x(2); cos(t.*w)-x(1).^3.*0.2-x(1)-x(2).*0.2)];
and in the first loop, you change from iterating over ‘e’ to iterating over ‘w’, with the values for ‘w’ as stated in the problem. I will leave those changes to you, since they are straightforward. The second loop changes similarly to reflect the changes in the problem requirements.
Sam
2015-7-12
Star Strider
2015-7-12
My pleasure!
I’m not quite sure what you mean by ‘green function’. (I am using R2015a, and it uses the parula colormap as its default. The lines were blue and red when I plotted them.) In my code, I plotted both the function ‘x(:,1)’ and the first derivative, ‘x(:,2)’.
To plot only the function and not the derivative, just plot ‘x(:,1)’, or equivalently (in my latest code that also detects the first zero-crossing), the ‘xv’ variable. Since I saved them as cells, it will probably be easier to create and plot ‘xv’.
saddam mollah
2018-7-9
编辑：saddam mollah
2018-7-9
Can you plot this problem in 3dimension as: t along x-axis, ev along y-axis and x1 along z-axis? Thank you in advance.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
