what the reason behind roundoff error inside while loop

2 次查看（过去 30 天）
rohit
rohit 2024-7-23
回答： Shubham 2024-7-23
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
  1 个评论
Aquatris
Aquatris 2024-7-23
编辑：Aquatris 2024-7-23
Ran in:
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

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Shubham
Shubham 2024-7-23
Ran in:
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.
Alan Stevens
Alan Stevens 2024-7-23
Ran in:
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

类别

在 Help Center 和 File Exchange 中查找有关 Multirate Signal Processing 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by