How do I use an If-statement when using RK4?

Below is a RK4 program for which I solve VM.

Here's my problem - in my assignment, there is a certain condition for which the values for VM(2) will be reduced by 70% and VM(3) by 20% each and every year, i.e VM(2)=VM(2)*0.3 and VM(3)=VM(3)*0.8 which my code below does not include since I do not know how to include it.

I do not know how to enter this into to the RK4 while loop and then plot the graf from t to t_end. I am assuming that I need to use an if-statement but not sure how to do it. Does anyone know how to enter my condition in my RK4-program and then how to plot the new values from t:t_end?

Grateful for any help that I can get!

a1=16;
a2=1.8*10^-5;
a3=0.011;
b1=2.0;
b2=0.085;
%Calculating my starting values for VM
v_start=(b1.^(10/6)*a1/a3)/(a2/a3*b1.^(10/6)+b2.^(10/6));
s_start=a1/a3-a2/a3*v_slut;

%T= nr of years
t_end=12; t=3; 
VM=[v_start, s_start 2]; T=t; VS=VM;
h=1/365;

while t<t_end-h/2
    f1=fvs2(t,VM);
    f2=fvs2(t+h/2, VM+h*f1/2);
    f3=fvs2(t+h/2, VM+h*f2/2);
    f4=fvs2(t+h, VM+h*f3);
    VM=VM+h/6*(f1+2*f2+2*f3+f4); t=t+h; T=[T; t]; VS=[VS;VM];
end,

subplot(3,1,1), plot(T, VS(:,1)) 
subplot(3,1,2),plot(T,VS(:,2))
subplot(3,1,3),plot(T,VS(:,3))

Answer 1

[Abraham Boayue]

You actually did not mention the condition. Can you say exactly what the condition is? We know what will happen, but we need to know the condition as well.

Answer 2

[Abraham Boayue]

Is this what you want to do? I am having a bit of a difficulty understanding what you want to do. I set v_sutt = 1, since there was no value assigned to it.

V1 = 0.3*VS(:,2);
V2 = 0.8*VS(:,3);
subplot(3,1,1)
plot(T,VS(:,1),'linewidth',2,'color','b');
a = title('Something Else');
set(a,'fontsize',14);
a = xlabel('T[ 3 12]');
set(a,'fontsize',14);
a = ylabel('values'); 
set(a,'fontsize',14);
grid
subplot(3,1,2)
plot(T,V1,'linewidth',2,'color','r');
a = title('30% of Pests population');
set(a,'fontsize',14);
a = xlabel('T[ 3 12]');
set(a,'fontsize',14);
a = ylabel('values'); 
set(a,'fontsize',14);
grid
subplot(3,1,3)
plot(T,V2,'linewidth',2,'color','g');
a = title('80% of Predictor population');
set(a,'fontsize',14);
a = xlabel('T[ 3 12]');
set(a,'fontsize',14);
a = ylabel('values'); 
set(a,'fontsize',14);
grid

[Tanmoy Bari]

Please find my answer below.

Answer 3

[Tanmoy Bari]

Sure.

The problem statement I need help solving follows below (translated from Swedish). Note that this is a continuation of a previous problem which I already solved with RK4. The problem is now how to enter the following statement into my code:

"At the time T3, they decide to spray poison so that 70 percent of the pests are killed at each year's spraying, which occurs at the end of the year. The effect is unfortunately that even the predators are affected, 20 percent of the predators are killed at the same time each year by poison."

Pests in my code is VM(2) and the Predators are VM(3). Therefore, I formulated the statement which I want to use as the following, as mentioned above:

VM(2)=VM(2)*0.3 VM(3)=VM(3)*0.8

Below is my function fvs2 if it is to any help.

      function res=fvs2(t,VM);
      V=VM(1); S=VM(2); R=VM(3);
      a1=16;
      a2=1.8*10^-5;
      a3=0.011;
      b1=2.0;
      b2=0.085;
      b3=1.5;
      c1=2;
      c2=0.025;

      res=[a1*V-a2*V^2-a3*V*S
          -b1*S^1.4+b2*V^0.6*S^0.8-1.5*S*R
          -c1*R+c2*S*R^0.5]';

Answer 4

[Tanmoy Bari]

Abraham,

Thank you so much for your reply and effort. I am sorry that I am not making the problem clear. But I'll try again. What I am aiming to do is to reduce the values: VM(2)=VM(2)*0.3 VM(3)=VM(3)*0.8 for each year, i.e. for t=3,4,5....12 and see how the growth looks like by plotting the graphs. So I want to insert this into my RK4-while loop.

The correct results should look like this:

Again, truly appreciate your help.

Answer 5

[Abraham Boayue]

Hi Tanmoy, sorry that I haven't been able to get back to you since our last communication.

In the following line of your code,

%Calculating my starting values for VM
v_start=(b1.^(10/6)*a1/a3)/(a2/a3*b1.^(10/6)+b2.^(10/6));
s_start=a1/a3-a2/a3*v_slut;

What is the value for v_slut?