Executing the while loop

Guys there is something wrong with my while loop. The while loop is supposed to run for 7389 iterations , but it only does two. what can be the problem ? Im just trying to compare the runge kutta method and the euler implicit method for a given ODE.

clear all;close all; clc;
format long;
u=1.5;
S=-33;
L=12/100;
N=11;
delx=L/(N-1);
C(1)=0.414;
x(1)=0;
%4th Order RK Method
for i=1:N-1
    x(i+1)=x(i)+delx;
    k1 = delx*(S/u*C(i));
    k2 = delx*(S/u*C(i)+S/u*k1/2);
    k3 = delx*(S/u*C(i)+S/u*k2/2);
    k4 = delx*(S/u*C(i)+S/u*k3);
    C(i+1) = C(i) + (1/6)*(k1+2*k2+2*k3+k4);
end
rk4=C(i+1);
ct_1=1;

% Euler Implicit
N_euler=2;
K(ct_1)=0.414;
y(ct_1)=0;
err=1;
while err>=0.00001
    delxeuler=L/(N_euler-1);
    y(ct_1+1)=y(ct_1)+delxeuler;
    K(ct_1+1)=(K(ct_1)/delxeuler)*(1/((1/delxeuler)-(S/u)));
    e(ct_1) = (rk4 - K(ct_1+1))^2;
    ct_1=ct_1+1;
    N_euler=N_euler+1;
    err = (sqrt(sum(e)))/N_euler;
end
plot(y,K,'-.+','color','g')
title('Concentration of CO vs. Distance');
xlabel('Axial(x) Direction [m]');
ylabel('Concentration of CO[mol/m3]');
hold on
plot(x,C,'-o','color','r')
legend('Euler Implicit: N=7389','Runge Kutta 4th Order: N=11');

[Image Analyst]

Did you know, in the MATLAB editor, you can type control-a (to select all) then type control-i to fix the indenting? It's hard to follow improperly indented code.

[DIP]

done!

Answer 1

[DIP]

Figured it out ! Thanks for your help ! 

% RK4 Method vs Implicit Euler Method

clear all;close all; clc;
format long;
u=1.5;
S=-33;
L=12/100;

%4th Order RK Method
N=11;
delx=L/(N-1);
C(1)=0.414;
x(1)=0;
for i=1:N-1
    x(i+1)=x(i)+delx;
    k1 = delx*(S/u*C(i));
    k2 = delx*(S/u*C(i)+S/u*k1/2);
    k3 = delx*(S/u*C(i)+S/u*k2/2);
    k4 = delx*(S/u*C(i)+S/u*k3);
    C(i+1) = C(i) + (1/6)*(k1+2*k2+2*k3+k4);
end
rk4=C(i+1);

% Euler Implicit Method

N_euler=1;
K(1)=0.414;
y(1)=0;
imp=0;

while abs(rk4-imp)>=0.00001
    N_euler = N_euler+1;
    delxeuler=L/(N_euler-1);

      for j=1:N_euler-1
          %Euler Implicit
          y(j+1)=y(j)+delxeuler;
          K(j+1)=(K(j)/delxeuler)*(1/((1/delxeuler)-(S/u)));
          imp=K(j+1);
          err=rk4-imp;
      end

end

plot(y,K,'-.+','color','g')
title('Concentration of CO vs. Distance');
xlabel('Axial(x) Direction [m]');
ylabel('Concentration of CO[mol/m3]');
hold on
plot(x,C,'-o','color','r')
legend('Euler Implicit: N=7389','Runge Kutta 4th Order: N=11');

[John BG]

cool, would you please be so kind to elaborate, just a little bit, and explain to the interested readers of your question the following?

1.

where was the error?

2.

why is it you wanted the amount of loops to reach precisely none other of the mentioned amount of 7389?

3.

how do you know in advance the amount of loops that RK or Euler will take to go below certain threshold, if a true test should be done with more or less random inputs, thus the error should vary, and then the amount of loops required to reach target should vary too?

thanks in advance for time and attention, awaiting any comment

John BG

<mailto:[e-mail] [e-mail]>

[Rahul Kalampattel]

Good job, glad you figured it out!

Answer 2

[Rahul Kalampattel]

When I run your code, after the first iteration err=0.028062920372021. Since this doesn't satisfy the condition err<=0.00001, the loop breaks.

I'm guessing the 'greater than' sign should be a 'less than' sign, since you want to iterate until convergence.

err=1; 

while err>=0.0001
   % do calcs, find error
end

[DIP]

sir, i made that change, but the calculation should end at 7389 iterations and not go on beyond that.

[Rahul Kalampattel]

There might be something wrong with your equations, or you might need to adjust your max error.

How do you know there should only be that number of iterations? If you want exactly 7398 iterations, I guess you could code it into the loop.

while N_euler<7390
 % do calcs, find error
end

[DIP]

[DIP]

i should be getting this.

[Rahul Kalampattel]

Try this, I took the spacing out of the loop.

err=1;
N_euler_max=7389;
delxeuler=L/(N_euler_max-1);

while N_euler<=N_euler_max
    y(ct_1+1)=y(ct_1)+delxeuler;
    K(ct_1+1)=(K(ct_1)/delxeuler)*(1/((1/delxeuler)-(S/u)));
    e(ct_1) = (rk4 - K(ct_1+1))^2;
    ct_1=ct_1+1;
    N_euler=N_euler+1;
    err = (sqrt(sum(e)))/N_euler;
end

(Edit: plot seems ok?)

[DIP]

Hi Rahul, I have to code the program in such a way that I am unaware of the final iteration. N=7389 is the answer im supposed to get.

[John BG]

Rahul

DIP wants the while loop stop at 7389 without the brake.

Indeed a loop may stop at a given amount of iterations if you limit such amount.

But DIP means that the error has to dive below 0.00001 at precisely 7389 laps.

[Rahul Kalampattel]

Sorry, but I'm not really sure what else to do. The way the code is at the moment, at 7389 iterations the error is greater than 0.00001 so the loop isn't going to stop. I guess you can try modifying the way you calculate your error term, changing one of these two lines:

e(ct_1) = (rk4 - K(ct_1+1))^2;
err = (sqrt(sum(e)))/N_euler;

Without knowing the actual equations you're trying to implement, I can't say for sure if there is a mistake somewhere else.

One more thing, unless the step size for the implicit Euler method is small enough, you're not going to a smooth plot. In your code right now, in the first iteration delxeuler=L/(N_euler-1)=L/(2-1)=L. If you plot this, you'll see it's nothing like the figure above. You either need to edit this, or take it out the loop and give it a constant value like I did here.

Answer 3

[John BG]

Hi DIP

I agree with Rahul

You can get the loop to spin more times by increasing the error threshold.

1.

..
err=1
while err>0.001
..

ct_1 stops at 872

tic toc measures 0.018 seconds

2.

    err=1
    while err>0.0001
    ..

ct_1 stops at 87305

tic toc just before and after the while loop measures 5.488 seconds

3.

err=1
while err>0.00001
..

I stopped it around a minute

ct_1 reached 631381

but the error still 3.7186e-5

it seems as if your coding of RK does not reduce the error below 0.00001 within a reasonable time, if it does it at all.

If you find this answer useful would you please be so kind to mark my answer as Accepted Answer?

To any other reader, please if you find this answer of any help solving your question,

please click on the thumbs-up vote link,

thanks in advance

John BG

<mailto:[e-mail] [e-mail]>

[DIP]

Sir, how can i improve the code of RK ? It is straightforward.