0

I tried loops of for and while with if inside but its showing only zeros in the n matrix , some how the loop is not getting incorporated. Please help me finding what is the error in this case particularly. also there is no error in the forumla and all please help me find the error in the computational process only. The error is probably in the for and nested for loop i hope

clear all;
clc
D = 1e-4 * 1e-4 ;                         % 1e-4 cm^2/sec
J=1e6 * 1e6;                              %input flux at half
format long;
%% Grid defination
del_x = 1e-9;
del_t=1e-6;
GridPoints_x=0:del_x:10e-6;
GridPoints_t=0:del_t:1e-5;
npoints_x=length(GridPoints_x);
npoints_t=length(GridPoints_t);
n=zeros(npoints_x,npoints_t);
Source = zeros(npoints_x, npoints_t);% Source Profile

flux_location_x = 5 * 1e-6;                  % Injection Location
flux_location_t = 1 * 1e-6;                  % Injection Time
f_index_x = find(GridPoints_x == flux_location_x);% Location Index
f_index_t = find(GridPoints_t == flux_location_t);% Time Index
i = 1;
Source(f_index_x, f_index_t) = J ;
tol = 1e-10;
dn=10;
while 1
n_old=n;
for m=2:npoints_t
    
    for i= 2:npoints_x-1
            n(i , m) = ((1/del_x^2) * (n(i + 1, m) + n(i - 1, m)) +n(i, m - 1)/(D * del_t) + Source(i, m) ...
    )/(2/del_x^2 + 1 / (D * del_t));
    end 
    n_new=n;
    if(abs(n_old-n_new)<=tol)
        break
    end
 end
end
Subham
  • 13
  • 4
  • I'm fairly certain your error is related to numerical precision, see [Why is 24.000 not equal to 24.000](https://stackoverflow.com/questions/686439/why-is-24-0000-not-equal-to-24-0000-in-matlab). Your `find_idx_x` is empty, whereas if you do `f_index_x = find(abs(GridPoints_x - flux_location_x)<1e-10)` it comes out as 5001, which is what you want. Since it's empty in your code, Source remains empty and the algorithm fails. Summary: don't naively compare numbers using `==`. – Adriaan Sep 05 '22 at 05:33
  • Also, you're making max. 1e5 iterations in your nested `for` loops every `while` run, changing only a single number the first time. This sounds inefficient. More efficient would be to e.g. check which of the `n` and `Source` entries are non-zero and only calculate that point, i.e. only calculate `i=5001; m=2` in the first iteration, rather than looping over loads and loads of zeros every time. – Adriaan Sep 05 '22 at 05:37
  • Final point: doing `while` loops with only a single condition, in your case a break upon reaching the tolerance, might get you stuck in an infinite loop (if the tolerance is never reached). I'd always include an iteration counter for `while` loops and also stop the loop on some predefined maximum iteration number, i.e. `while iter < maxiter; iter = iter + 1; end` to make sure it stops on `maxiter` iterations, as well as your tolerance condition. – Adriaan Sep 05 '22 at 06:20
  • Thank you everyone I got the point of == operatpr and did the necessary changes in the while loop as well. I did not know the problem with == operator. – Subham Sep 05 '22 at 09:53
  • One minor detail: you can pre-calculate several factors in your `for` loops. `(1/del_x^2)`, `(D * del_t)` and especially `(2/del_x^2 + 1 / (D * del_t))` can be pre-calculated before the loops, given they are constant. They're basic operations, but doing those millions of times will take some additional time. – Adriaan Sep 05 '22 at 13:44

0 Answers0