#C19.txt, April 3, 2014 #Numerical solutions of the Heat Equation Help:=proc(): print(`Heq1(U,x,h,x1,t1), Heq1f(U,x,h,x1,t1)`): print(` Heq2(U,x,h,x1,t1) `): end: #Heq1(U,x,h,k,x1,t1): Explicit method #inputs an expression U of the variable x #two small numbers h and k (the mesh sizes in x and t resp.) #and a final time t1 #outputs an APPROXIMATION to the value of u(x1,t1) #with x1 and t1 multiples of h and k. #u(x,t) is the solution of the 1D-Heat equation #u_{xx}=u_t, u(x,0)=U(x), u(0,t)=0, u(1,t)=0 Heq1:=proc(U,x,h,k,x1,t1) local T,j,n,N,J,r: N:=t1/k: J:=1/h: r:=k/h^2: for j from 0 to J do T[j,0]:=subs(x=j*h,U): od: for n from 0 to N-1 do T[0,n+1]:=0: T[J,n+1]:=0: for j from 1 to J-1 do T[j,n+1]:=(1-2*r)*T[j,n]+r*T[j-1,n]+r*T[j+1,n]: od: od: #T[i,j] is an approximation to u(i*h,k*j) T[x1/h,t1/k]: end: #Heq1f(U,x,h,k,x1,t1): floating version of Heq1(U,x,h,k,x1,t1): Heq1f:=proc(U,x,h,k,x1,t1) local T,j,n,N,J,r: N:=t1/k: J:=1/h: r:=k/h^2: for j from 0 to J do T[j,0]:=evalf(subs(x=j*h,U)): od: for n from 0 to N-1 do T[0,n+1]:=0: T[J,n+1]:=0: for j from 1 to J-1 do T[j,n+1]:=evalf((1-2*r)*T[j,n]+r*T[j-1,n]+r*T[j+1,n]): od: od: #T[i,j] is an approximation to u(i*h,k*j) T[x1/h,t1/k]: end: #Heq2(U,x,h,k,x1,t1): Implicit method same input and output as in Heq1(U,x,h,k,x1,t1): #inputs an expression U of the variable x #two small numbers h and k (the mesh sizes in x and t resp.) #and a final time t1 #outputs a list of lists such that #L[i][j] is an APPROXIMATION to the value of u(x1,t1) #of the temperature at x=i*h and t=j*k #u(x,t) is the solution of the 1D-Heat equation #u_{xx}=u_t, u(x,0)=U(x), u(0,t)=0, u(1,t)=0 Heq2:=proc(U,x,h,k,x1,t1) local T,j,n,N,J,r,eq,var: N:=t1/k: J:=1/h: r:=k/h^2: for j from 0 to J do T[j,0]:=subs(x=j*h,U): od: for n from 0 to N-1 do var:={seq(T[j,n+1], j=0..J)}: eq:={ T[0,n+1]=0, T[J,n+1]=0}: eq:= eq union {seq((1+2*r)*T[j,n+1]-r*T[j-1,n+1]-r*T[j+1,n+1]=T[j,n],j=1..J-1)}: var:=solve(eq,var): for j from 0 to J do T[j,n+1]:=subs(var,T[j,n+1]): od: od: #T[i,j] is an approximation to u(i*h,k*j) T[x1/h,t1/k]: end: