Why "solve" gives a different result, if everything generates the same?

15 views Asked by At

First code:

X := 8*sqrt(x);
u0 := 0;
s0 := 1/2;
x0 := 3/5;
N := 300;
h := 10/N;
E0 := 3/5;
E1 := 3.5;
R := proc (i) options operator, arrow; E0*Heaviside(i+x0-x)+E1*Heaviside(x-i-x0) end proc;
i -> E0 Heaviside(i + x0 - x) + E1 Heaviside(x - i - x0)
for i from 0 to N do
  xi[i] := h*i 
end do;
for j to N-1 do 
     b[j] := evalf((Int((x-xi[j-1])*X, x = xi[j-1] .. xi[j])+Int((xi[j+1]-x)*X, x = xi[j] .. xi[j+1]))/h) 
end do;
b[N] := evalf((Int((x-xi[N-1])*X, x = xi[N-1] .. xi[N]))/h);
for i to N do
  for m to 10 do
   if m-1 <= xi[i] and xi[i] <= m then 
        E[i] := evalf((int(R(m-1), x = xi[i-1] .. xi[i]))/h) 
   end if 
  end do 
end do;

eq[0] := u[0]-u0;
eq[N] := evalf((-E[N-1]*u[N]+E[N-1]*u[N-1])/h+b[N]+s0);
for i to N-1 do eq[i] := evalf((E[i]*u[i-1]-(E[i]+E[i+1])*u[i]+E[i+1]*u[i+1])/h+b[i]) end do;
EQ := {seq(eq[i], i = 0 .. N)};
U := {seq(u[i], i = 0 .. N)};
A := solve(EQ, U);

Second code:

trapezoidalRule := proc (E, a, b, n, k) 
local h, result, x, i;
h := (b-a)/n; 
result := (1/2)*E(a, k-1)+(1/2)*E(b, k-1); 
for i to n-1 do 
   x := a+i*h; 
   result := result+E(x, k-1) 
end do; 
result := result*h; 
return result 
end proc;
for i from 0 to N do 
  xi[i] := h*i 
end do;
for j to N-1 do 
  b[j] := evalf((Int((x-xi[j-1])*X, x = xi[j-1] .. xi[j])+Int((xi[j+1]-x)*X, x = xi[j] .. xi[j+1]))/h) 
end do;
b[N] := evalf((Int((x-xi[N-1])*X, x = xi[N-1] .. xi[N]))/h);
for i from 0 to N do
    if i = 0 then 
         eq[i] := u[i]-u0 
    elif i = N then 
       _local(lastE); 
       for j to 10 do 
          if E(xi[i], j-1) <> 0 then
             lastE := evalf(trapezoidalRule(E, xi[i-2], xi[i-1], N, j)/h) 
          end if
       end do; 
       eq[i] := evalf((-lastE*u[N]+lastE*u[N-1])/h+b[N]+s0) 
     else 
        _local(firstE, secondE);
        for j to 10 do 
           if E(xi[i], j-1) <> 0 then
            firstE := evalf(trapezoidalRule(E, xi[i-1], xi[i], N, j)/h)
           end if; 
           if E(xi[i+1], j-1) <> 0 then
            secondE := evalf(trapezoidalRule(E, xi[i], xi[i+1], N, j)/h) 
           end if 
        end do; 
        eq[i] := evalf((firstE*u[i-1]-(firstE+secondE)*u[i]+secondE*u[i+1])/h+b[i])
      end if 
end do;
EQ := {seq(eq[i], i = 0 .. N)};
U := {seq(u[i], i = 0 .. N)};
A := solve(EQ, U);

And for some reason this "resolve" it gives different results for them, although these variables in it are absolutely the same. I will be glad if someone helps to figure it out.

I checked all variables, they are the same for 2 cases.enter image description hereenter image description here[enter image description here] enter image description here(https://i.stack.imgur.com/oehAL.png)

0

There are 0 answers