Reputation: 5
restart;
N := 100;
h := 10/N;
xi := Vector(N+1);
for i to N+1 do
xi[i] := (i-1)*h
end do;
X := proc (x) options operator, arrow; 8*sqrt(x) end proc;
Xi := Vector(N+1);
for i to N+1 do
Xi[i] := evalf(X(xi[i]))
end do;
s0 := 1/2;
E := proc (x, i) options operator, arrow;
piecewise(i <= x and x <= i+3/5, 3/5, i+3/5 < x and x <= i+1, 3.5, 0)
end proc;
E0 := E(xi[1], 0);
u := Vector(N+1);
u0 := 0;
for i to N+1 do for j to 10 do
if i = 1 then
u[i] := 0
elif `and`(i = N+1, E(xi[i-1], j-1) <> 0) then
EN := (E(xi[i-1], j-1)+E(xi[i-2], j-1))*(1/2);
u[i] := (-E0*v[i]+EN*v[i-1])/h+(1/2)*h*Xi[i]+s0
else
firstEi; secondEi;
if E(xi[i-1], j-1) <> 0 then
firstEi := (E(xi[i], j-1)+E(xi[i-1], j-1))*(1/2)
end if;
if E(xi[i], j-1) <> 0 then
secondEi := (E(xi[i+1], j-1)+E(xi[i], j-1))*(1/2)
end if;
u[i] := (firstEi*v[i-1]-(firstEi+secondEi)*v[i]+secondEi*v[i+1])/h+h*Xi[i]
end if
end do;
A := Matrix(N+1, N+1);
for i to N+1 do
for j to N+1 do
if `and`(i = 1, j = 1) then
A[i, j] := 0 else A[i, j] := coeff(u[i], v[j])
end if
end do
end do;
B := Vector(N+1);
for i to N+1 do
B[i] := coeffs(u[i])[1]
end do;
triangleMatrix := gaussjord(concat(A, B));
solution := backsub(triangleMatrix);
R := [seq([xi[n], solution_reverse[n]], n = 1 .. N+1)];
listplot(R, color = red);
I implemented the variation-difference method and wanted to display a graph, but I don't understand where the error is if the graph should be as I drew it in green, and I get this:
Well, that is, the graph is correct, but it must be turned around, if I am not mistaken, the sign of the second derivative affects it.
Upvotes: 0
Views: 44