what the reason behind roundoff error inside while loop

while t<=1

I2=exp(-pi*(t-0.5*dt))*I1;

Z=(inv(Bu))*F1;

B1v=2*A/dt+B-F-F2*Z;

B12v=(F+F2*Z-B+2*A/dt)*v0 -2*I2;

v1=B1v\B12v;

q=F1*(v0+v1);

q1=Bu\q;

u1=(q1-u0);

u1=u1;

v0=v1;

u0=u1;

% x=linspace(l,L,n);x1=x(2:n-1);

uex_val= uex(x1, t)';

vex_val = vex(x, t)';

t=t+dt;

[t,dt]

end and dt=(0.1)*(16^(-p)) ,for p=1 ' when i am running this loop over t i should get at last t=1 , but somehow after so many steps this toop is taking an approximate value of dt =.00624999 instead of .00625 what could be the possible reason for that

[Aquatris]

I do not see see the problem. I removed the unncessary parts which do not change t or dt and run your code and got expected results. You might be have defined t and dt as global variables and changing them within the uex() and/or vex() functions, if they are functions.

t = 0;
p = 1;
dt=(0.1)*(16^(-p));
while t<=1
      t=t+dt;
end
format long
t

t = 
   1.006249999999997

dt

dt = 
   0.006250000000000

Answer 1

[Shubham]

Hi Rohit,

The issue you're encountering is likely due to the accumulation of floating-point errors when incrementing t by dt in each iteration of the loop. Floating-point arithmetic can introduce small errors because not all decimal fractions can be represented exactly in binary. Over many iterations, these small errors can accumulate, leading to a final value of t that is slightly off from the expected value.

% Parameters
p = 1; % Example value for p
dt = (0.1) * (16^(-p));
t_final = 1;
num_steps = round(t_final / dt);

% Initial conditions
t = 0;
n = 10; % Example size of the problem (adjust as needed)
v0 = ones(n, 1); % Example initial value for v0 (adjust as needed)
u0 = ones(n, 1); % Example initial value for u0 (adjust as needed)

% Define other necessary variables and functions
I1 = ones(n, 1); % Example value for I1 (adjust as needed)
Bu = eye(n); % Example value for Bu (adjust as needed)
F1 = ones(n, 1); % Example value for F1 (adjust as needed)
A = eye(n); % Example value for A (adjust as needed)
B = eye(n); % Example value for B (adjust as needed)
F = eye(n); % Example value for F (adjust as needed)
F2 = eye(n); % Example value for F2 (adjust as needed)
x1 = linspace(0, 1, n); % Example x1 (adjust as needed)
x = linspace(0, 1, n); % Example x (adjust as needed)

% Define example functions uex and vex
uex = @(x, t) sin(pi * x) * exp(-t); % Example function (adjust as needed)
vex = @(x, t) cos(pi * x) * exp(-t); % Example function (adjust as needed)

for step = 1:num_steps
    I2 = exp(-pi * (t - 0.5 * dt)) * I1;
    Z = inv(Bu) * F1;
    B1v = 2 * A / dt + B - F - F2 * Z;
    B12v = (F + F2 * Z - B + 2 * A / dt) * v0 - 2 * I2;
    v1 = B1v \ B12v;
    q = F1 .* (v0 + v1); % Element-wise multiplication
    q1 = Bu \ q;
    u1 = q1 - u0;
    
    % Update variables
    v0 = v1;
    u0 = u1;
    
    % Calculate exact values
    uex_val = uex(x1, t)';
    vex_val = vex(x, t)';

    % Update time
    t = t + dt;
end

% Display final time and time step
disp([t, dt]);

    1.0000    0.0063

The loop runs for a fixed number of steps (num_steps), which is calculated based on the total time (t_final) and the time step (dt). This approach avoids the issue of floating-point error accumulation in the time variable t.

Notes:

• Adjust the size n and any initial values or functions to match your specific problem requirements.
• Ensure that any matrix operations are dimensionally consistent with the data you are using.

Answer 2

[Alan Stevens]

Compare the following

dt = 0.1*16^-1;

t = 0;
while t<1
    t = t+dt;
end
disp(t)

    1.0062

t = 0;
c = 1;
while t<1
    t = c*dt;
    c = c+1;
end
disp(t)

     1