Question: Need to reduce the Jacobian matrix construction as number of equations increases

Dear Users!
Hope everyone fine here. I tried (as given bellow) to find the solution of nonlinear system of PDEs via FDM. To solve system of nonlinear equations I used newton raphson method and for higher value of like Mx > 8 the matrix G and G1 (mentioned as red) take alot of time. Can anyone help me to reduce the computational time? Becuase I have to evealuate the solution for Mx = 50.


restart; Digits := 30; with(LinearAlgebra);
T := 1; L := 3; N := 30; Mx := 5; `Δx` := L/(1.*Mx); `Δt` := T/(1.*N);
for i from 0 while i <= Mx do
u[i, 0] := 0.; u[i, -1] := u[i, 1]; tau[i, 0] := 0.; theta[i, 0] := 0.; theta[i, -1] := theta[i, 1]
end do;
for n from 0 while n <= N do u[0, n] := 0.; u[Mx, n] := 0.; theta[0, n] := 1.; theta[Mx, n] := 0.
end do;
for n from 0 while n <= N-1 do
print("Simulation in proccess at time-level n", n+1);
for i while i <= Mx-1 do
Ru[i, n] := simplify((u[i+1, n+1]-u[i+1, n])/`&Delta;t`+(u[i+1, n+1]-2*u[i+1, n]+u[i+1, n-1])/`&Delta;t`^2-(u[i+1, n+1]-2*u[i, n+1]+u[i-1, n+1])/`&Delta;x`^2+25.*(u[i+1, n+1]+(u[i+1, n+1]-u[i+1, n])/`&Delta;t`)-1.5*(theta[i, n]+(theta[i, n+1]-theta[i, n])/`&Delta;t`));
`R&theta;`[i, n] := simplify((theta[i+1, n+1]-theta[i+1, n])/`&Delta;t`+(theta[i+1, n+1]-2*theta[i+1, n]+theta[i+1, n-1])/`&Delta;t`^2-(theta[i+1, n+1]-2*theta[i, n+1]+theta[i-1, n+1])/((15.)*`&Delta;x`^2)-((u[i, n+1]-u[i-1, n+1])/`&Delta;x`)^2/(3.)) end do;
for i while i <= Mx-1 do
`R&tau;`[i, n] := simplify(tau[i+1, n+1]+(tau[i+1, n+1]-tau[i+1, n])/`&Delta;t`-1.5^(-1/4)*(u[i+1, n+1]-u[i, n+1])/`&Delta;x`)
end do;
Sys := `<,>`(seq(Ru[i, n], i = 1 .. Mx-1), seq(`R&tau;`[i, n], i = 1 .. Mx-1), seq(`R&theta;`[i, n], i = 1 .. Mx-1));
V := `<,>`(seq(u[i, n+1], i = 1 .. Mx-1), seq(theta[i, n+1], i = 1 .. Mx-1), seq(tau[i, n+1], i = 2 .. Mx));
G := Matrix(3*(Mx-1), proc (i, j) options operator, arrow; diff(Sys[i], V[j]) end proc); G1 := MatrixInverse(G);
X[n, 0] := Vector(1 .. 3*(Mx-1), 1);
for k1 from 0 to r do
X[n, k1+1] := eval(V-G1 . Sys, Equate(V, X[n, k1]))
end do;
Sol[n] := Equate(V, X[n, r+1]); assign(op(Sol[n]));
if n > 0 then
U := eval(`<,>`(seq(u[i1, n+1], i1 = 1 .. Mx)-seq(u[i1, n], i1 = 1 .. Mx))); Noru[n+1] := Norm(%, 2); print("L[&infin;] norm of &tau;(x,y,t) at time level = ", %);
Theta := eval(`<,>`(seq(theta[i1, n+1], i1 = 0 .. Mx)-seq(theta[i1, n], i1 = 0 .. Mx))); `Nor&theta;`[n+1] := Norm(%, 2); print("L[&infin;] norm of &theta;(x,y,t) at time level = ", %) else print("n < 0")
end if end do

Special request to:
@acer @Carl Love @Kitonum @Preben Alsholm

Please Wait...