Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Replies: 1   Last Post: Jan 14, 2013 11:27 PM

 Messages: [ Previous | Next ]
 Nasser Abbasi Posts: 6,677 Registered: 2/7/05
Posted: Jan 14, 2013 11:27 PM

On 01/13/2013 11:01 PM, Fred Bartoli wrote:
> Hello,
>
> As a first simple test case of a more complicated pb I want to solve the
> plane diffusion equation.
> It seems NDSolve gives a bad answer to the simple step case but I can't
> find what's happening...
>
> (* Set up the equation: *)
>
> eq1=D[h[t,x],x,x]- k1 D[h[t,x],t]==0
>
> const=k1-> 10000
>
> (*"solve" it...*)
>
> sol = NDSolve[
> {D[h[t, x], {x, 2}] - k1 D[h[t, x], t] == 0, h[0, x] == 0,
> h[t, 1] == 0, h[t, 0] == 1 - Exp[-10^6*t]} /. const,
> h, {t, 0, 10^-3}, {x, 0, 0.01}][[1]]
>
> Plot3D[(h[t,x])/.mag/.solCart,{x,0,10^-4},{t,0,10^-4},AxesLabel->{"x","t","h"},PlotRange->All,PlotLabel->"Bogus
> solution"]
> Plot3D[Erfc[x/(2*Sqrt[t/k1])]/.const,{x,0,10^-4},{t,0,10^-4},AxesLabel->{"X","t","h"},PlotRange->All,PlotLabel->"The
> right one"]
>
> (* Test the Erfc[x/(2Sqrt[t/k1])] solution *)
> Release[Hold[D[h[t,x],x,x]- k1*D[h[t,x],t]]/.
> h[t,x]->Erfc[x/(2*Sqrt[t/k1])]]//Simplify
>
>
> I hope I'm doing something wrong, but I can find what...
>

First, your code, little cleaned up to make it easy to
see and to remove the `mag` and `solCart` which I have
no idea what it is. There is no need to write `eq=`
but later on, copy the same equation again. Why?

Also, it is best to write the IC in one variable,
and the equations in another, so things are easier
to see.

----------------------------------------------
eq1 = D[h[t, x], x, x] - k1 D[h[t, x], t] == 0
ic = {h[0, x] == 0, h[t, 1] == 0, h[t, 0] == 1 - Exp[-10^6*t]}
const = k1 -> 10000;

sol = First@
NDSolve[Flatten@{eq1, ic} /. const, h, {t, 0, 1}, {x, 0, 0.01}];

plot[f_, x_, t_, title_] :=
Plot3D[f, {x, 0, 0.01}, {t, 0, 1}, AxesLabel -> {x, t, h},
PlotRange -> All, PlotLabel -> title]

plot[h[t, x] /. sol, x, t, "Bogus solution"]
plot[Erfc[x/(2*Sqrt[t/k1])] /. const, x, t, "The right one"]
------------------------------------------

The question though, why do you think that the initial
conditions will not make any difference to the solution
of the pde? You are comparing a numerical solution for some
arbitrary IC you specficed above with one particular solution
that you have found.

--Nasser

Date Subject Author
1/14/13 Fred Bartoli
1/14/13 Nasser Abbasi