Help in writing a function

Hi everyone, I need to write a function to solve this system:

Until now all I can came up with was:

 function dz = control1(v,z,parameters)

% gammaT=parameters(1);
phi_0=parameters(2);
psi_0=parameters(3);
psi_c0=parameters(4);
B=parameters(5);
Lc=parameters(6);
W=parameters(7);
H=parameters(8);
C = 0;

gammaT_max = 0.9;
gammaT_min = 0.61;

A = (gammaT_max - gammaT_min)/2; %amplitude
b = gammaT_max - ((gammaT_max - gammaT_min)/2);
hertz = 50;
w=hertz*2*pi;
gammaT = @(t) A*sin(w*t)+b 

dz = zeros(2,1);
dz(1)=(1/(4*B*B*Lc))*(z(2)-gammaT(v)*(z(1))^0.5);
psi_c=psi_c0+H*(1+1.5*(z(2)/W-1)-0.5*(z(2)/W-1).^3);
dz(2) =(1/Lc)*(psi_c-z(1));
 end
 

which is obviously wrong since gammaT = @(t) A*sin(w*t)+b shoouldn't be define this way.

I think gammaT should be written in that way but the time in the argument of the sin should be rearranged in some way I don't know

Answer 1

[Star Strider]

‘which is obviously wrong since gammaT = @(t) A*sin(w*t)+b shoouldn't be define this way.’

Why? It appears to be coded correctly with respect to the posted image. It’s being evaluated with respect to your independent variable ‘v’, that appears to be correct. If it should actually be a different variable (such as time), it would be necessary to define the time in the context of the existing variables. We would need more information in order to help you with that.

For what it’s worth, the function runs without error for me using:

parameters = rand(8,1);

and:

[V,Z] = ode45(@(v,z)control1(v,z,parameters), [0 10], rand(2,1));

figure
plot(V, Z)
grid

.

[Paul Rogers]

exactley, it's like yoou said, I need to be like

gammaT = (A*sin(w*t)+b)

where t is the time.

I attached al the files I use too run this function in the previous post

[Star Strider]

Before I go through all the files, do you have a PDF of the system you’re modeling that you can attach here? (Reading it would likely be easier for me to understand, since I seriously doubt that what you’re doing is anywhere close to my areas of expertise.)

I’d like to see how ‘v’ and ‘t’ are related, or if they even are.

[Paul Rogers]

yes, the sistem is the one described in eq. 1 and 2 but in eq. 1 Phi_t(Psi) = gamma_T*(Psi)^0.5.

In my case gammaT is a function of the time described by the equation:

gammaT = (A*sin(w*t)+b)

I forgot to link the article, it's just the system in the equation 1 and 2 withe the modification about gammaT

Greitzer System

I can't upload the files I am using because I've reached the limit of 10 files, but they can be found in my previous post here:

https://www.mathworks.com/matlabcentral/answers/621798-problem-with-a-function#comment_1077308

[Star Strider]

Thank you.

It will take me a bit to review that to understand how time relates to everything else. Compressors and the equations modeling them are entirely new to me.

EDIT — (31 Oct 2020 at 15:19)

Equation (6) may be important here:

since ξ is non-dimensional time. The equations for:

and:

would imply (to me) that they are functions indirectly of time, so whatever you are using to represent time (is it ‘v’?) would appear to me to be the argument to the function. .

I need a bit of clarification.

[Paul Rogers]

I thank you, it's really important to me.

Anyway I did the functioons in the paper, but now too add something from me I need this small adjstment about gammaT, wich is noot a constant anymore but a function of the time

gammaT = A*sin(w*t)+b

[Star Strider]

My poleasure!

Please see my previous Comment after EDIT that I just posted.

[Star Strider]

I cannot find ‘gammaT’ in the paper, at least not by that name.

Can you reference an equation where it appears?

[Paul Rogers]

gammaT in the paper is in table 1 and its value is 0.61.

The problem is I don't have to use gammaT as a constant like in the paper, but in my case gammaT varies with a form llike:

gammaT = A*sin(w*t)+b

where:

gammaT_max = 0.9;

gammaT_min = 0.61;

A = (gammaT_max - gammaT_min)/2; %amplitude

b = gammaT_max - ((gammaT_max - gammaT_min)/2);

hertz = 50;

w=hertz*2*pi;

[Star Strider]

Are the values (range) of ‘A’ and ‘b’ between the limits (‘gammaT_min’ and ‘gammaT_max’) functions of or controlled by anything else, or any other variables or parameters (since ‘w’ appears to be fixed)?

[Paul Rogers]

no, it's like you said, A’ and ‘b’ are between the limits (‘gammaT_min’ and ‘gammaT_max’),

[Star Strider]

I’m a bit lost, here. The ‘gammaT’ function appears to be well-defined. In the code you originally posted, the code for ‘control1’ runs without error and appears to give acceptable results.

Is ‘gammaT’ supposed to be something else — some other function — or behave differently from what you coded it to be? Iy is certainly possible to define it as a constant if you want it to be constant, or to be something different from a sine curve (perhaps atan or tanh)?

I must be missing something, because I still don’t understand what the problem is.

[Paul Rogers]

yes, I tought it was correct as I defind gammaT in control1, and the result looked acceptabe to me too.

The problem is in reality it's wrong, (my proofessor told me). When I plot this part (from general)

figure()
plot(time,psi,'b','Linewidth',3);

xlabel('Time (s)')
ylabel('$\Psi$','Interpreter','latex')
grid

I expect to see a function who has fluctuations like a sin and where I could see the 50Hz, but all I gget is a line with no meaning.

If in line 32 and 33 from general, I swithc from ode45 to ode113

[time1 ans1] = ode113(@(v,z)nocontrol(v,z,compressor_parameters),[0 t_c/timeRatio],[0 0]);
[time2 ans2] = ode113(@(v,z)control1(v,z,compressor_parameters),[t_c/timeRatio t_f/timeRatio],[real(ans1(end,1)) real(ans1(end,2))]);
       

I get a line which has the shape I expected (like a sinusoid), but still don't see the 50Hz. All I can count is 3 or 4 peaks in 2 seconds which is roughly 1.5-2.0Hz.

[Paul Rogers]

so, with ode113 I get the shape I expected but stil I don't see the 50Hz.

At this point, I am starting to think that since I need to "see" 50Hz the time-step I need is 10ms at least, but when I open the time vector I se that time-step is at least one order more.

Now the problem could be this, but I really have no clue in how to improve the time discretization of an ode

[Paul Rogers]

I think I got it, so, the function and stuff are well written, the prooblem is now the time step that the ode45 and the oders use. it's too big, infact if I simulate a low frequency, i.e. 0.001, I get the curve I expect.

Now I think I should open another topic asking how to improve the time resolution in an oode

[Walter Roberson]

make your tspan a vector of 3 or more values instead of 2 values. The internal time step will be the same but you can can get results to whatever resolution you ask for.

[Paul Rogers]

thanks, but how can I do this, I am really lost at this point.

this is where I call the function:

[time1 ans1] = ode113(@(v,z)nocontrol(v,z,compressor_parameters),[0 t_c/timeRatio],[0 0]);
[time2 ans2] = ode113(@(v,z)control1(v,z,compressor_parameters),[t_c/timeRatio t_f/timeRatio],[real(ans1(end,1)) real(ans1(end,2))]);
    

I really have no idea where to start now.

Here another call of the same function (it has a different name but the ssystem is the same):

[t,y]=ode113(@greitzer_Jerzak,[0 1200],[0,0]);

[Star Strider]

If you want the time to go from 0 to 1200 with a step of 0.001:

tspan = linspace(0, 1200, 1200000);

will work. (I prefer linspace to the colon operator because it’s more accurate.)

[Paul Rogers]

how do I modify this:

t,y]=ode45(@greitzer_Jerzak,[0 40],[0,0]);

[Star Strider]

The [0 40] ‘tspan’ argument would be:

tspan = linspace(0, 40 , 40000);

to return solutions in steps of 0.001. The third argument here is simply the second argument x1000, to get that specific resolution.

[Paul Rogers]

it's incredible, I don't understand.

If I put a frequency up to 0.1Hz I get the correct result, when I go to higher frequencies (50Hz) the system it's not working anymore.

[Star Strider]

I don’t understand how that could be the situation, other than that you may be seeing ‘aliasing’.

One other option is to increase the length of the ‘tspan’ vectors, for example multiplying the third argument by 1000 to decrease the step size (using linspace) by 1000.

[Paul Rogers]

if I put

hertz = 0.1; 

in the function and

tspan = linspace(0, 100 , 600000);
[t,y]=ode113(@greitzer_Jerzak,[tspan],[0,0]);

in the main, I exacley get the result I expected, but if I want to investigate a higher frequency more appropriate for my studies, like 50hz, it doesn't work

[Walter Roberson]

tspan = linspace(0, 40 , 40000);

should be

tspan = linspace(0, 40 , 40000+1);

because the endpoints are included in the count.

[Paul Rogers]

tspan = linspace(0, 100 , 600000);
[t,y]=ode113(@greitzer_Jerzak,[tspan],[0,0]);

with this I have an incredible resultion, like 0.0001 in the time step, but it's stilll shoowing aliasing at 50Hz.

[Paul Rogers]

thank you @Walter Roberson, but it doesn't change a thing. always good up to 0.1hz, but noot right for 50Hz

[Paul Rogers]

the problem with the solution:

tspan = linspace(0, 100 , 600000);
[t,y]=ode113(@greitzer_Jerzak,[tspan],[0,0]);

is that the function is evalueted on the fixed points I specified, but the algoritm is still choosing its own points to solve the system.

[Star Strider]

That’s how the MATLAB ODE solvers work.

They're adaptive internally with respect to the steps, and a vector of 3 or more elements for ‘tspan’ reports the outputs at the specified values of the independent variable.The ‘tspan’ vector otherwise doesn’t affect the internal workings of the ODE solvers, at least as I understand them.

[Paul Rogers]

sure, but is there a way to solve my problem at 50Hz?

[Paul Rogers]

here is the solution I was looking for, so I can observe up tp 50Hz:

init=[0 0]'; 
options= odeset('MaxStep',0.001); %maximum time-step size
[t,y]=ode45(@greitzer_Jerzak,[0,20],init,options);

This solution allows me to chose the maximum step size of the time, evven if ode45 will still pick up its own points, but at least I make sure they are at a sample frequeny much higher so I can catch the frequencies I need.

thank you everybody, expecially to Star Strider who really put a lot of efford to point me in the right direction.

[Star Strider]

As always, my pleasure!

(I can Comment again!)