# Solve analytically 1D transport equation?

Posted 5 months ago
948 Views
|
11 Replies
|
3 Total Likes
|
 I'm trying to solve for the analytical solution of 1D transport equation to verify the results of the numerical solution. eqn = D[u[x, t], t] == D[u[x, t], {x, 2}]; ic = u[x, 0] == 50; bc1 = u[0, t] == 50; bc2 = D[u[1, t], x] == 0; DSolve[{eqn, ic, bc1, bc2}, u[x, t], {x, t}] However, I obtain the following, DSolve::deqn: Equation or list of equations expected instead of True in the first argument {(u^(0,1))[x,t]==(u^(2,0))[x,t],u[x,0]==50,u[0,t]==50,True}. Have I missed any step? I'm looking for the symbolic solution of the PDE with Dirichlet boundary condition at the inlet and Neumann boundary condition at the outlet. Could someone help?
11 Replies
Sort By:
Posted 5 months ago
 You need to write the no flux BC differently, but your specification returns a trivial result of 50 when corrected. eqn = D[u[x, t], t] == D[u[x, t], {x, 2}]; ic = u[x, 0] == 50; bc1 = u[0, t] == 50; bc2 = (D[u[x, t], x] == 0) /. {x -> 1}; DSolve[{eqn, ic, bc1, bc2}, u[x, t], {x, t}] (* {{u[x, t] -> 50}} *) 
Posted 5 months ago
 Many thanks for the response.I would like to if it is possible to view the complete expression of the analytical solution. I think the solution that is displayed after making the above-suggested change is the steady-state solution.
Posted 5 months ago
 I think that you will need to rethink your specification. Right now, your initial condition says that the temperature is 50 everywhere. Your first boundary condition says that the temperature at the boundary is 50. Your second boundary conditions says that the system is insulated at the end. How can the system change from the value of 50 everywhere?
Posted 5 months ago
 Solving for analytical solution of 1D transport equation
Posted 5 months ago
 Solving for analytical solution of 1D transport equation
Posted 5 months ago
 Thanks for the advice. I have changed the left boundary condition bc1 = u[0, t] == 70; sol = DSolveValue[{eqn, ic, bc1, bc2}, u[x, t], {x, t}] gives, 70 + 2 Inactive[Sum][( 40 E^(-225 \[Pi]^2 t (-1 + 2 K[1])^2) Sin[1/2 \[Pi] x (-1 + 2 K[1])])/(\[Pi] - 2 \[Pi] K[1]), {K[1], 1, \[Infinity]}] Since the above solution is an infinite series, I am considering first 300 terms. asol = sol /. {\[Infinity] -> 300} Next, I am trying to plot the solution at x = 0.5 for time t= 1:20 using the following command Plot(sol) The above syntax isn't complete, I would like to ask for some help in specifying the syntax for plotting the above solution.
Posted 5 months ago
 Look at the DSolveValue>Applications>Classical Partial Differential Equations for an example how to do this. eqn = D[u[x, t], t] == D[u[x, t], {x, 2}]; ic = u[x, 0] == 50; bc1 = u[0, t] == 70; bc2 = (D[u[x, t], x] == 0) /. {x -> 1}; dsol = DSolveValue[{eqn, ic, bc1, bc2}, u[x, t], {x, t}] /. {K[1] -> m}; asol[x_, t_] = dsol /. {\[Infinity] -> 300} // Activate; Plot[asol[0.5, t], {t, 1, 20}, PlotRange -> All] Note that a value of 300 causes a warning about the value being too small.
Posted 5 months ago
 Thank you. I made the following changes, eqn = D[u[x, t], t] == 100*D[u[x, t], {x, 2}]; and asol[x_, t_] = dsol /. {\[Infinity] -> 30000} // Activate; Plot[asol[0.5, t], {t, 1, 20}, PlotRange -> All] However, the plot shows only a straight line at u = 70 (y-axis) from 0 to 20 seconds.I had expected the plot to start from 50 at t = 0 and reach 70 after the initial transient phase. I'm not sure what's causing the problem. Any suggestions on how to get the right solution?
Posted 5 months ago
 Hi Natasha,Your diffusion coefficient is huge and that implies that the system responds very quickly. To see this response, you will need to look at times much closer to 0 and not begin at $t=1$. You probably do not need many terms to capture the response. eqn = D[u[x, t], t] == 100 D[u[x, t], {x, 2}]; ic = u[x, 0] == 50; bc1 = u[0, t] == 70; bc2 = (D[u[x, t], x] == 0) /. {x -> 1}; dsol = DSolveValue[{eqn, ic, bc1, bc2}, u[x, t], {x, t}] /. {K[1] -> m}; asol5[x_, t_] = dsol /. {\[Infinity] -> 5} // Activate; asol10[x_, t_] = dsol /. {\[Infinity] -> 300} // Activate; Plot[{asol5[0.5, t], asol10[0.5, t]}, {t, 0, 0.001}, PlotRange -> All] 
 Thank you very much. I varied the value diffusion constant and check the response.I want to check how the analytical solution differs when there is convection along with diffusion. eqn = D[u[x, t], t] == 100*D[u[x, t], {x, 2}] - 50*D[u[x, t], x] ; For the same set of initial and boundary conditions, I'm checking for the analytical solution dsol = DSolveValue[{eqn, ic, bc1, bc2}, u[x, t], {x, t}] /. {K[1] -> m}; The output that I obtain is , DSolveValue[{ \!$$\*SuperscriptBox[\(u$$, TagBox[ RowBox[{"(", RowBox[{"0", ",", "1"}], ")"}], Derivative], MultilineFunction->None]\)[x, t] == -100 \!$$\*SuperscriptBox[\(u$$, TagBox[ RowBox[{"(", RowBox[{"1", ",", "0"}], ")"}], Derivative], MultilineFunction->None]\)[x, t] + 50 \!$$\*SuperscriptBox[\(u$$, TagBox[ RowBox[{"(", RowBox[{"2", ",", "0"}], ")"}], Derivative], MultilineFunction->None]\)[x, t], u[x, 0] == 50, u[0, t] == 70, \!$$\*SuperscriptBox[\(u$$, TagBox[ RowBox[{"(", RowBox[{"1", ",", "0"}], ")"}], Derivative], MultilineFunction->None]\)[1, t] == 0}, u[x, t], {x, t}] The summation terms are missing. Could you please help?
 Hi Natasha,You are probably approaching the limit of DSolve's ability to find an analytical solution (at least directly). You should consider a numerical solution for this problem.If you are okay with a numerical solution, you can use the code below as a starting point. tmax = 2/100; ptmax = tmax/5; eqn = D[u[x, t], t] == 100*D[u[x, t], {x, 2}] - 50*D[u[x, t], x]; ic = u[x, 0] == 50; bc1 = u[0, t] == 50 + 20 Tanh[20000 t] ; bc2 = (D[u[x, t], x] == 0) /. {x -> 1}; ndsol = NDSolveValue[{eqn, ic, bc1, bc2}, u[x, t], {x, 0, 1}, {t, 0, tmax}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 50}}, MaxStepFraction -> 1/1000]; ContourPlot[ndsol, {t, 0., ptmax}, {x, 0., 1}, PlotPoints -> {200, 200}, MaxRecursion -> 4, PlotTheme -> "Web", ColorFunction -> ColorData["DarkBands"]] Note that I used a Tanh function to provide a rapidly rising "step" function so that the boundary and initial conditions were consistent. If you want to explore the effects of velocity on the solution, you could wrap the above (note that flipped the plot 90 degrees so it would flow left to right) into a Module like so. condiff[v_] := Module[{tmax, ptmax, eqn, ic, bc1, bc2, ndsol, cp}, tmax = 2/100; ptmax = 2 tmax/5; eqn = D[u[x, t], t] == 100*D[u[x, t], {x, 2}] - v*D[u[x, t], x]; ic = u[x, 0] == 50; bc1 = u[0, t] == 50 + 20 Tanh[20000 t] ; bc2 = (D[u[x, t], x] == 0) /. {x -> 1}; ndsol = NDSolveValue[{eqn, ic, bc1, bc2}, u, {x, 0, 1}, {t, 0, tmax}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 2000}}, MaxStepFraction -> 1/1000]; cp = ContourPlot[ndsol[x, t], {x, 0., 1}, {t, 0., ptmax}, PlotRange -> {50, 70}, PlotPoints -> {50, 50}, MaxRecursion -> 4, PlotTheme -> "Web", ColorFunction -> ColorData["DarkBands"], FrameLabel -> Automatic]; Show[{cp, Graphics[Arrow[{{0, 0.004}, {v/500, 0.004}}]]}] ] Now, you can create an animation of velocities from 0 to 500 (it will take a while!). imgs = Table[condiff[v], {v, 0, 500, 5}]; ListAnimate[imgs]