How to update variable according to solutions?

Hi, I am almost new in matlab. I am trying to update my variable according to error calculation. I wrote my code by using ıf statement but I cannot found how can I update my data if my error is not in the range of my tolerance. Also my code is not working now. I assigned an error but I guess my if statement is not true :( Can anyone tell me how can I update my dt after error calculation? This is my code,

t0=0;
h0=100;
tEnd = 1000;
dt =50;
N=(tEnd-t0)/dt;

T = [t0:dt:tEnd];
Y = zeros(N+1,1);
Y(1) = h0;
%loop using euler's method
for i = 1:N
    fi = 20 - sqrt(Y(i));
    Y(i+1) = Y(i) + fi * dt ; 
end
eu = Y(i+1); %answer by using euler
plot(T,Y)
hold on
% loop for modified euler's method
for i = 1:N
    k1 = 20 - sqrt(Y(i));
    tNew = T(i) + dt;
    yNew = Y(i) + dt*k1;
    k2 = 20 - sqrt(yNew);
    Y(i+1) = Y(i) + (0.5 * dt * (k1+k2)) ; 
end
plot(T,Y)
meu = Y(i+1); %answer by using modified euler
e = 10^-2 ; %error
if abs(meu-eu)*0.5*dt <=e
    disp('error is in the range of tolerance');
else
    return 
    dt

end
fprintf ( ' %f, %f \n' , eu,meu);

[Geoff Hayes]

PÜRLEN - please clarify what you mean by your code is not working now. I can run it and get an answer with two solutions (from the Euler and Modified Euler algorithms) that are roughly the same but don't meet your tolerance requirement. Is that the problem? Or is the problem that you have a return call in your else statement

if abs(meu-eu)*0.5*dt <=e
    disp('error is in the range of tolerance');
else
    return 
    dt
end

and so the final line (the fprintf) is not being called?

[PÜRLEN SEZER]

The problem is the error is not below my error tolerance. It calculated euler and modified euler parts but I cannot return dt value and update the value. Since I don't know how to change dt value after finding that the error is higher than my tolerance, I said it is not working.

[Geoff Hayes]

How would you change dt if the tolerance is not met? Would you increase or decrease this value?

If you want to the code to continue running with a new dt, then you could put the majority of your (above) code into a while loop that is conditioned on the tolerance check

 meu = [];
 eu = [];
 dt = 50;
 e = 10^-2 ;
 while isempty(meu) || abs(meu-eu)*0.5*dt > e 
    % update dt
    dt = ???;
    % rest of your code
 end 
 disp('error is in the range of tolerance');
 fprintf ( ' %f, %f \n' , eu,meu);

[PÜRLEN SEZER]

I want to decrease my dt. You mean that I have to write while loop for dt value at the beginning of my code and rest of the code remains same. But I don't know what isempty(Meu) does? I am really new at matlab. Sorry, if my questions are not clear for you?

[Roger Stafford]

@PÜRLEN: If Geoff’s ‘while’ coding is not clear to you, an equivalent form could be this:

 dt = 50;
 e = 10^(-2);
 b = true;
 while b
   % your code computes eu and meu
   b = abs(meu-eu) >= e;
   if b
     % decrease dt (better still, increase N)
   end 
 end

You will have to decide how you want to decrease dt. Perhaps divide it by 2? It will need to be such that (tEnd-t0)/dt remains an integer. It might be easier to simply increase N and compute dt = (tEnd-t0)/N.

[It is clear you are working with the differential equation:

dy/dt = 20-sqrt(y), y(0) = 100

In case you are interested, this does have a solution in terms of the implicit equation:

t = -2*sqrt(y)-40*log(20-sqrt(y))+20+40*log(10)

which can be expressed explicitly with y a function of t using the lambertw function.]

[PÜRLEN SEZER]

Thanks for your answers but I cannot still change my dt value after error calculation is done. I want to change it after it is found bigger than my tolerance value. And this is final version of my code.

t0=0;
h0=100;
tEnd = 1000;
dt =50;
N=(tEnd-t0)/dt;

T = [t0:dt:tEnd];
Y = zeros(N+1,1);
Y(1) = h0;
%loop using euler's method

for i = 1:N
    fi = 20 - sqrt(Y(i));
    Y(i+1) = Y(i) + fi * dt ; 
end
eu = Y(i+1); %answer by using euler
plot(T,Y)
hold on
% loop for modified euler's method
for i = 1:N
    k1 = 20 - sqrt(Y(i));
    tNew = T(i) + dt;
    yNew = Y(i) + dt*k1;
    k2 = 20 - sqrt(yNew);
    Y(i+1) = Y(i) + (0.5 * dt * (k1+k2)) ; 
end

meu = Y(i+1)  %answer by using modified euler   

fprintf ( ' %f, %f \n' , eu,meu);

e = 10^-2 ; %error tolerance

while abs(meu-eu)*0.5*dt > e

t0=0;
h0=100;
tEnd = 1000;
dt =1 ; %update dt
N=(tEnd-t0)/dt;
T = [t0:dt:tEnd];
Y = zeros(N+1,1);
Y(1) = h0;

      for i = 1:N
      fi = 20 - sqrt(Y(i));
      Y(i+1) = Y(i) + fi * dt ; 
  end
  eu = Y(i+1); %answer by using euler
  plot(T,Y)
  hold on
  % loop for modified euler's method
  for i = 1:N
      k1 = 20 - sqrt(Y(i));
      tNew = T(i) + dt;
      yNew = Y(i) + dt*k1;
      k2 = 20 - sqrt(yNew);
      Y(i+1) = Y(i) + (0.5 * dt * (k1+k2)) ; 
  end

end

Answer 1

[Geoff Hayes]

PÜRLEN - your code is fixing dt to be one on each iteration of the while loop

while abs(meu-eu)*0.5*dt > e
  t0=0;
  h0=100;
  tEnd = 1000;
  dt =1 ; %update dt

You need to update dt. For example, your code could become

atIteration = 0;
while abs(meu-eu)*0.5*dt > e
  atIteration = atIteration + 1;
  t0=0;
  h0=100;
  tEnd = 1000;
  if atIteration > 1
      dt = dt + 1; % or do whatever you want to increase dt
  end

The above will change your dt but the change might not have a positive effect. And so as you increase dt, you may diverge from the best found solution. I suspect that you will need some code to handle that scenario.

[Geoff Hayes]

Your code now has some duplication. You shouldn't need to redo the Euler calculation on each iteration of the while loop and you can do the modified Euler calculation only within the while loop (no need to do this outside as well).