Advanced Search

Question: Improve my code for parabolic equations

Question:Improve my code for parabolic equations

mailman 81 Maple
plot format
October 25 2007
0
I have been told my code is not 'perfect', I'm not too hot on maple and it has taken me a week to get to this stage. Any help or suggestions would be very welcome! h:=0.1: k:=0.001: r:=k/h^2: xl:=0: xr:=1: # left and right ends of interval in x nxpoints:=round((xr-xl)/h): ntsteps:=5: # total number of time steps iprint:=1: # print output every iprint time steps printf(cat(" x,t ","%7.3f"$(nxpoints+1),"\n"),seq(i*h,i=0..nxpoints)): format:=cat("%7.3f"$(nxpoints+2),"\n"): # format for printing output Boundary conditions for j to ntsteps+1 do u[0,j]:=0: u[nxpoints,j]:=0 end do: Initial condition for i from 0 to nxpoints do u[i,0]:=sin(Pi*i*h) end do: solution:=[[0.,seq(u[i,0],i=0..nxpoints)]]: printf(format,op(solution[1])): Time steps for j from 0 to ntsteps-1 do for i from 1 to nxpoints-1 do u[i,j+1]:=r*u[i-1,j]+(1-2*r)*u[i,j]+r*u[i+1,j] end do; solution:=[op(solution),[(j+1)*k,seq(u[i,j+1],i=0..nxpoints)]]: if floor((j+1)/iprint)=(j+1)/iprint then printf(format,op(solution[j+2])) end if; end do: Plots for j from 1 by iprint to ntsteps+1 do plot||j:=plot([seq([i*h,solution[j,i+2]],i=0..nxpoints)],labels=[x,u]): end do: plots[display]({seq(plot||(iprint*j+1),j=0..ntsteps/iprint)}); Comparison with exact solution plottime:=0.005: uexact:=(x,t)->sin(Pi*x)*exp(-Pi^2*t): pexact:=plot(uexact(x,plottime),x=0..xr,labels=[x,u],linestyle=DASH): plots[display]({pexact,plot||(floor(plottime/k)+1)});
Please Wait...