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])/`Δt`+(u[i+1, n+1]-2*u[i+1, n]+u[i+1, n-1])/`Δt`^2-(u[i+1, n+1]-2*u[i, n+1]+u[i-1, n+1])/`Δx`^2+25.*(u[i+1, n+1]+(u[i+1, n+1]-u[i+1, n])/`Δt`)-1.5*(theta[i, n]+(theta[i, n+1]-theta[i, n])/`Δt`));

`Rθ`[i, n] := simplify((theta[i+1, n+1]-theta[i+1, n])/`Δt`+(theta[i+1, n+1]-2*theta[i+1, n]+theta[i+1, n-1])/`Δt`^2-(theta[i+1, n+1]-2*theta[i, n+1]+theta[i-1, n+1])/((15.)*`Δx`^2)-((u[i, n+1]-u[i-1, n+1])/`Δx`)^2/(3.)) end do;

for i while i <= Mx-1 do

`Rτ`[i, n] := simplify(tau[i+1, n+1]+(tau[i+1, n+1]-tau[i+1, n])/`Δt`-1.5^(-1/4)*(u[i+1, n+1]-u[i, n+1])/`Δx`)

end do;

Sys := `<,>`(seq(Ru[i, n], i = 1 .. Mx-1), seq(`Rτ`[i, n], i = 1 .. Mx-1), seq(`Rθ`[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[∞] norm of τ(x,y,t) at time level = ", %);

Theta := eval(`<,>`(seq(theta[i1, n+1], i1 = 0 .. Mx)-seq(theta[i1, n], i1 = 0 .. Mx))); `Norθ`[n+1] := Norm(%, 2); print("L[∞] norm of θ(x,y,t) at time level = ", %) else print("n < 0")

end if end do

Special request to:

@acer @Carl Love @Kitonum @Preben Alsholm