what the reason behind roundoff error inside while loop
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaire
Aquatris
le 23 Juil 2024
Modifié(e) : Aquatris
le 23 Juil 2024
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
dt
Réponses (2)
Shubham
le 23 Juil 2024
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]);
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.
0 commentaires
Alan Stevens
le 23 Juil 2024
Compare the following
dt = 0.1*16^-1;
t = 0;
while t<1
t = t+dt;
end
disp(t)
t = 0;
c = 1;
while t<1
t = c*dt;
c = c+1;
end
disp(t)
0 commentaires
Voir également
Catégories
En savoir plus sur Multirate Signal Processing dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!