Asked by dcydhb dcydhb
on 11 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

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.

dcydhb dcydhb
on 12 Jun 2019

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

dcydhb dcydhb
on 11 Jun 2019

Can you be more specific?thanks！

Rik
on 11 Jun 2019

See my editted answer.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.