Variable is in a function file, how do I use a loop to change the value of this variable?
3 次查看(过去 30 天)
显示 更早的评论
If the variable is in the main program, I can compile a loop to change the value of a variable and get the result with the variable, but what if the variable is in the function file?
it means i want to use loop to get the plot that zmax changes as the cd or cm changes.
main program are as this
clc;
clear all;
close all;
step=0.1;
t0=50;
length1=t0/step;
opts = odeset('RelTol',1e-4,'AbsTol',1e-4);
Z=zeros(length1,2);
[T,Z]=ode45('vdp', [0:step:t0], [0,1], opts);
figure(1)
plot(T,Z(:,1),'r.');
zend=length(Z(:,1));
zmax=max(Z(20:zend,1))%get the plot that zmax changes as the cd or cm changes
function file is as this
function dz = vdp(t,z)
g=9.8;
gamma=10094;
d=20;
J0=15;
J1=5;
c1=7680;
r=0.21;
m=13500;
G=m*g;
c2=32000;
m1=13500;
m2=0;
%************************************************************
%************************************************************
cd=0.8;%which I want to use the loop to change
%************************************************************
%************************************************************
k0=2000;
p0=1000;
area=15.9;
%************************************************************
%************************************************************
ca=1;%which I want to use the loop to change
%************************************************************
%************************************************************
h=0.5;
t0=3.5;
d1=0.8283479061205706;
height=1;
rho=1025;
dz=zeros(2,1);
dz(1)=z(2);
if z(2)>=0
a1x=1;
else a1x=0;
end
if z(2)>=0
Jx=J0;
else Jx=J1;
end
if z(2)>=0
b1x=0;
else b1x=0;
end
dz(2)=...
...
-(p0/3 + g*m + (k0*z(1))/9 - area*g*rho*(d1 - z(1) + (h*cos((2*pi*t)/t0))/2) + (a1x*c2*z(2))/r + (area*cd*rho*abs(z(2) + (pi*h*sin((2*pi*t)/t0))/t0)*(z(2) + (pi*h*sin((2*pi*t)/t0))/t0))/2 + (2*pi^2*area*h*rho*cos((2*pi*t)/t0)*(d1 - z(1) + (h*cos((2*pi*t)/t0))/2))/t0^2 + (2*pi^2*area*ca*h*rho*cos((2*pi*t)/t0)*(d1 - z(1) + (h*cos((2*pi*t)/t0))/2))/t0^2)/(m + Jx/r^2 + area*ca*rho*(d1 - z(1) + (h*cos((2*pi*t)/t0))/2));
end
0 个评论
回答(2 个)
Jan
2019-6-11
Do not redefine the important Matlab function cd as a variable, because this can cause serious troubles during debugging.
[T,Z]=ode45('vdp', [0:step:t0], [0,1], opts);
Providing a function to be integrated as a char vector is outdated for over 15 years now. Although this is still supported for backwar compatibility, it is less powerful and secure than a function handle. In addition a function handle allows to define a parameter:
p = 3.1415; % Your parameter
[T,Z] = ode45(@(t,y) vdp(t, y, p), 0:step:t0, [0,1], opts);
function dy = vdp(t, y, p)
...
end
By the way, o:step:to is a vector already. Using the Matlab operator for array concatenation [] is a waste of time here only.
Another problem can cause invalid outputs: Your function to be integrated is not smooth:
if z(2)>=0
Jx=J0;
else Jx=J1;
end
Matlab's ODE integrators are designed to process smooth functions only. Otherwise the step size controller can drive completely mad, which can cause different effects: A stop with an error message (minimal step size exceeded), or even worse a result, which is dominated by rounding errors due to a massive increase of the evaluated steps. Don't do this. Do not use a numerical tool for a computation, which does not match the specified limitations.
Use event functions instead, which stop the integration when the state of the function changes.
Rik
2019-6-11
编辑:Rik
2019-6-11
According to the doc you can pass extra parameters by setting an anonymous function. I see no indication in the R2014a doc that this wouldn't work in that release.
If you make the changes indicated below, you should be able to have cd and ca as inputs, so you can change them in a loop.
cd=0.7;
ca=1.1;
[T,Z]=ode45(@(t,z) vdp(t,z,cd,ca), [0:step:t0], [0,1], opts);
function dz = vdp(t,z,cd,ca)
if nargin<4
%set defaults
ca=1;
if nargin<3
cd=0.8;
end
end
%%
%rest of your function
%%
end
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!