Executing a function in iterations

Problem Statement: I have a system of differential equations in which i need to run a particular function for n iterations using preferably for loop. So basically i want my main script to run for all the different values of function.

I have used a similar model equation in this question and just posting my attempted code here.

I must mention that i have MATLAB 2016a in which i can't use function file in the main script.

function:MWE_fn.m

function rk1 = MWE_fn(t,y)
r2=0.5;b2=1;d_A=0.05; 
 c=0.5;
rk1(1)=r2*y(1)*(1-b2*y(1))-c*y(1)*y(2);
rk1(2)= A0(t)- d_A*y(2);
rk1=rk1(:);
end

function fa=A0(t)
    if t>round(t)
        fa=0.4;
    else 
        fa=0;
    end
    
end

Script: MWE_body.m

timerange= 0:0.5:360;
IC= [0.1,0];%initial conditions 
for idx=1:length(timerange)-1
    [tTemp,yTemp] =ode45(@(t,y) MWE_fn(t,y),timerange(idx:idx+1), IC);
    t(idx) = tTemp(end); 
    y(idx,:) = yTemp(end,:);
    IC = y(idx,:);
end
plot(t,y(:,1),'b');
hold on
plot(t,y(:,2),'r');

The above is an working example and i attempted to use a for loop in the function file but it showed me the error that i cannot define a function like that. What i want is that the following function, should run for different values of fa like from 0.1 to 1 and give plots for all those values separately.

function fa=A0(t)
    if t>round(t)
        fa=0.4;
    else 
        fa=0;
    end
    
end

Thank you for any help you can provide.

Answer 1

[Walter Roberson]

Your formulas are discontinuous. You cannot use any of the ode* functions for them. The discontinuities occur exactly at the integer times and nowhere else.

If I recall correctly you asked about the same system with different phrasing before, talking about a system that is driven the first half of the day and not driven the second half. You were advised then on how to proceed, and you were advised then that ode* routines evaluate at intermediate times, not just the boundaries you give for tspan. For example ode45 will try to evaluate at roughly t=1.08 and your test is not prepared that; 1.08>round(1.08) yes but 1.0>round(1.0) false but 1.0 is part of the first half day and so should be forced.

Do not use an if in an ode routine. Instead have two different functions, one with forcing and one without, and alternate calling them for the half days.

[Vira Roy]

Thanks Walter for your reply.

Yes i did ask for the similar system few days back and you helped me there. I am not exactly very experienced in MATLAB and hence the trivial problems.

Can i not get away with the false statement you mentioned above 1.0>round(1.0) false by putting an equal to sign 1.0>=round(1.0) false?

Sorry for being blunt but do you mean i should use a different method of solving other than ODE* for this system? RK-4 or others would work? Also my original problem remains that, how do i loop the function?

Any help would be nice.

Thank you

[Walter Roberson]

timerange= 0:0.5:360;
IC= [0.1,0];%initial conditions 
for idx=1:2:length(timerange)-1
    [tTemp,yTemp] =ode45(@(t,y) MWE_fn_driven(t,y),timerange(idx:idx+1), IC);
    t(idx) = tTemp(end); 
    y(idx,:) = yTemp(end,:);
    IC = y(idx,:);
    [tTemp,yTemp] =ode45(@(t,y) MWE_fn_notdriven(t,y),timerange(idx:idx+1), IC);
    t(idx+1) = tTemp(end); 
    y(idx+1,:) = yTemp(end,:);
    IC = y(idx,:);
end
plot(t,y(:,1),'b');
hold on
plot(t,y(:,2),'r');

function rk1 = MWE_fn_driven(t,y)
r2=0.5;b2=1;d_A=0.05; 
 c=0.5;
rk1(1)=r2*y(1)*(1-b2*y(1))-c*y(1)*y(2);
rk1(2)= 0,4 - d_A*y(2);
rk1=rk1(:);
end
    
function rk1 = MWE_fn_notdriven(t,y)
r2=0.5;b2=1;d_A=0.05; 
 c=0.5;
rk1(1)=r2*y(1)*(1-b2*y(1))-c*y(1)*y(2);
rk1(2)= - d_A*y(2);
rk1=rk1(:);
end

[Vira Roy]

I thought my triviality led you away from the problem, but thank you for the working example. I don't want to push the help you gave, but if you can help with this final issue i would really be thankful.

That value of 0.4 you have used in the driven case, i want it to be a variable and need to see comapirison plot for different values, like 0.1, 0.2 etc.

i tried to modify the driven function as follows, but it didn't work

function rk1 = MWE_fn_driven(t,y)
for i=1:10
    A0(i)=i/10;
r2=0.5;b2=1;d_A=0.05; 
 c=0.5;
rk1(1)=r2*y(1)*(1-b2*y(1))-c*y(1)*y(2);
rk1(2)= A0(i) - d_A*y(2);
rk1=rk1(:);
end
end 
    

If you can show a code, it would be good else, even any kind of advise how to proceed would be fine.

P.S: I am using MATLAB 2016a, so i had to make separate function file for the functions in your modified code.

Thanks.

[Walter Roberson]

I am not at my desk at the moment so I cannot pull out the exact link, but search mathworks.com for "parameterize functions"

[Vira Roy]

Thanks Walter, the documentation on parametrization of functions really helped and i am accepting your answer.

Here is what i did.

function rk1 = MWE_fn_driven(t,y,A0)
r2=0.5;b2=1;d_A=0.05; 
 c=0.5;
% rk1(1)=r2*y(1)*(1-b2*y(1))-c*y(1)*y(2);
rk1(1)=r2*y(1)*(1-b2*y(1))-c*y(1)*y(2);
rk1(2)= A0- d_A*y(2); %Here i added the A0
 rk1=rk1(:);
end

main script

timerange= 0:0.5:360;
IC= [0.1,0.1];%initial conditions 
 for i=0.1:0.1:1 %Started the loop
     A0=i;
   for idx=1:2:length(timerange)-1
       
    [tTemp,yTemp] =ode45(@(t,y) MWE_fn_driven(t,y,A0),timerange(idx:idx+1), IC);
    t(idx) = tTemp(end); 
    y(idx,:) = yTemp(end,:);
    IC = y(idx,:);
    [tTemp,yTemp] =ode45(@(t,y) MWE_fn_notdriven(t,y),timerange(idx:idx+1), IC);
    t(idx+1) = tTemp(end); 
    y(idx+1,:) = yTemp(end,:);
    IC = y(idx,:);
    
end
plot(t,y(:,1),'b');
hold on
plot(t,y(:,2),'r');
IC= [0.1,0.1];%Reset the Initial conditions 
end 

Walter i just had two confusions

1. Is there a more elegant way to do this than this where i need to reset the Initial conditions ?
2. why does the graph thickness keeps increasing?

[Walter Roberson]

Is there a more elegant way to do this than this where i need to reset the Initial conditions ?

That depends on what you mean by "elegant". If you wanted, you could make your state vector be one (y and y') pair for each A0 value. There is no requirement that the state variables interact with each other. For example,

y = reshape(y, 2, :);
rk1(1,:)=r2*y(1,:)*(1-b2*y(1,:))-c*y(1,:)*y(2,:);
rk1(2,:)= A0 - d_A*y(2,:);
rk1=rk1(:);

for vector A0.

[Walter Roberson]

why does the graph thickness keeps increasing?

Your graph scale is too small to see the oscillation clearly. Remember you drive for 1/2 and do not drive for 1/2 so you expect a rapid oscillation.

[Vira Roy]

That clears everything Walter. Thanks a lot.