MATLAB Answers

0

Variable is in a function file, how do I use a loop to change the value of this variable?

Asked by dcydhb dcydhb on 11 Jun 2019
Latest activity Commented on by dcydhb dcydhb on 12 Jun 2019
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 Comments

Sign in to comment.

Products


Release

R2014a

2 Answers

Answer by Jan
on 11 Jun 2019

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.

  1 Comment

why we can't use the
if z(2)>=0
Jx=J0;
else Jx=J1;
end
and can you clarify the event functions usage by an example?
thanks a lot!!!

Sign in to comment.


Answer by Rik
on 11 Jun 2019
Edited by Rik
on 11 Jun 2019

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

  2 Comments

Sign in to comment.