清空筛选条件
清空筛选条件

what the reason behind roundoff error inside while loop

36 次查看（过去 30 天）
rohit
rohit 2024-7-23，8:25
回答： Shubham 2024-7-23，10:50
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，8:52
编辑：Aquatris 2024-7-23，8:53
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，10:50
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，8:53
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 中查找有关 Solver Outputs and Iterative Display 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by