# Executing the while loop

20 views (last 30 days)
DIP on 18 Feb 2017
Commented: Rahul Kalampattel on 20 Feb 2017
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');
##### 2 CommentsShowHide 1 older comment
DIP on 19 Feb 2017
done!

DIP on 20 Feb 2017
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');
##### 2 CommentsShowHide 1 older comment
Rahul Kalampattel on 20 Feb 2017
Good job, glad you figured it out!

Rahul Kalampattel on 19 Feb 2017
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
Rahul Kalampattel on 19 Feb 2017
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.

John BG on 19 Feb 2017
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.
John BG
DIP on 19 Feb 2017
Sir, how can i improve the code of RK ? It is straightforward.