regular pertubation for diffusion into a stagnant layer with reaction.

We can easily solve the case of a linear reaction expression (1st order)

DSolve[{D[θ[η],{η,2}] - ϵ θ[η]^1 == 0,θ[0] ==1, θ'[1] ==0},θ[η],η]

However, Mathematica crashes if the reaction is 2nd order

DSolve[{D[θ[η],{η,2}] - ϵ θ[η]^2 == 0,θ[0] ==1, θ'[1] ==0},θ[η],η]

However, Mathematica redeems itself with numerical solutions....

ans = NDSolve[{D[θ[η],{η,2}] - .01 θ[η]^2 == 0,θ[0] ==1, θ'[1] ==0},θ[η],η]

Plot[(θ[η]/.ans),{η,0,1},PlotLegends->{"ϵ=.01"},PlotStyle->Red]

ans = NDSolve[{D[θ[η],{η,2}] - .02 θ[η]^2 == 0,θ[0] ==1, θ'[1] ==0},θ[η],η]

Plot[(θ[η]/.ans),{η,0,1},PlotLegends->{"ϵ=.02"},PlotStyle->Green]

ans = NDSolve[{D[θ[η],{η,2}] - .05 θ[η]^2 == 0,θ[0] ==1, θ'[1] ==0},θ[η],η]

Plot[(θ[η]/.ans),{η,0,1},PlotLegends->{"ϵ=.05"},PlotStyle->Orange]

ans = NDSolve[{D[θ[η],{η,2}] - .2 θ[η]^2 == 0,θ[0] ==1, θ'[1] ==0},θ[η],η]

Plot[(θ[η]/.ans),{η,0,1},PlotLegends->{"ϵ=.2"},PlotStyle->Blue]

ans = NDSolve[{D[θ[η],{η,2}] - .5 θ[η]^2 == 0,θ[0] ==1, θ'[1] ==0},θ[η],η]

Plot[(θ[η]/.ans),{η,0,1},PlotLegends->{"ϵ=.5"},PlotStyle->Yellow]

numericalplot=Show[%45,%47,%49,%51,%53,PlotRange->All]

ans = NDSolve[{D[θ[η],{η,2}] - 100 θ[η]^2 == 0,θ[0] ==1, θ'[1] ==0},θ[η],η]

The numerical solution (that was done last week) also works for large “ϵ”, but you may recall that for very fast reactions, even the 1st order reaction caused some problems.

Plot[(θ[η]/.ans),{η,0,1},PlotLegends->{"ϵ=100"},PlotStyle->Dashing[{.01,.02}],PlotRange->All]

Here is the perturbation solution from Deen’s book.

perturbans = 1+ϵ/2( η^2-2 η) + ϵ^2 /3 (η^4/4 -η^3+ 2 η)

perturbplot =Plot[perturbans/.{ϵ->{.01,.02,.05,.2,.5}},{η,0,1}]

Compare with the numerical solutions...

Show[%,%109]

So there really is good agreement at sufficiently small ϵ which diminishes as ϵ gets bigger.  This is the expected result.

So the perturbation technique allowed us to get a usable analytical solution to a nonlinear ODE.  We could make it better by including more terms, but this usually reaches a point of limiting returns quickly.

perturbansO2 = 1+ϵ/2( η^2-2 η)

perturbplotO2 =Plot[perturbansO2/.{ϵ->{.01,.02,.05,.2,.5}},{η,0,1},PlotStyle->Dashing[{.05,.05}]]

We see that the O2 error version disagrees with the O3 error solutio after ϵ gets bigger than ~.05.   

We could do another higher term if we wished.

Show[perturbplotO2,perturbplot]

We could also compare the 2nd order kinetics to first order kinetics....

ans01 =DSolve[{D[θ[η],{η,2}] - ϵ θ[η]^1 == 0,θ[0] ==1, θ'[1] ==0},θ[η],η]

Plot[(θ[η]/.ans01)/.{ϵ->.2},{η,0,1},PlotLegends->{"ϵ=.2"},PlotStyle->Red]

Here is the numerical solution at ϵ = 0.2.

%106

Show[%,%106]

So the perturbation answer is probably better than assuming 1st order!

Can we solve the perturbation problem ourselves?
Here is the form of the solution

perturbexpansion = θ0[η] + ϵ θ1[η] + ϵ^2 θ2[η]+ ϵ^3 θ3[η]

The nonlinear ODE is

ode =  ϵ θ[η]^2 -     D[θ[η],{η,2}]

Here is how to substitute for the dependent variable and any derivatives

result1 =ode/.{θ[η] ->perturbexpansion,Derivative[a1_][θ][η]:>  (∂_({η,a1}))perturbexpansion}

We need to collect the separate powers of ϵ.  We need each power of ϵ to be separately equal to 0 otherwise, the results will be correct, or not correct depending on a specific value of ϵ.  We need an answer for all ϵ

result2 = Collect[ExpandAll[result1],ϵ]

order0 = Coefficient[result2,ϵ,0]

order1 =Coefficient[result2,ϵ,1]

order2 = Coefficient[result2,ϵ,2]

order3 = Coefficient[result2,ϵ,3]

Let’s get the anwers

order0answer =DSolve[{order0 ==0, θ0[0] == 1, θ0'[1] == 0}, θ0[η],η]

(θ0[η]^2-(θ1)''[η])/.order0answer[[1]]

order1

(%/. order0answer[[1]])==0

DSolve[1-(θ1)''[η]==0,{θ1[η]},{η}]

((order1 )/. order0answer[[1]] )==0

order1answer = DSolve[{((order1 )/. order0answer[[1]]) ==0, θ1[0] == 0, θ1'[1] == 0}, θ1[η],η]

order2answer = DSolve[{(((order2 )/. order1answer[[1]])/.order0answer[[1]]) ==0, θ2[0] == 0, θ2'[1] == 0}, θ2[η],η]

order3answer = DSolve[{((((order3 )/. order2answer[[1]])/.order1answer[[1]])/.order0answer[[1]]) ==0, θ3[0] == 0, θ3'[1] == 0}, θ3[η],η]

threetermsanswer = (((perturbexpansion/.(order3answer[[1]])/. order2answer[[1]])/.order1answer[[1]])/.order0answer[[1]])

perturbplotO2 =Plot[perturbansO2/.{ϵ->{.01,.02,.05,.2,.5}},{η,0,1},PlotStyle->Dashing[{.05,.05}]]

perturbploto3 =Plot[threetermsanswer/.{ϵ->{.01,.02,.05,.2,.5}},{η,0,1},PlotStyle->Dashing[{.01,.01}]]

Show[%,%%]

Compare the 3 term solution to the numerical result

numericalplot

Show[perturbploto3,numericalplot]

So with 1 more term we get to ϵ=.2 with no problem.  However, 0.5 is still not very good across the domain.  You could do another term if you wish!

Ok, one more term:

perturbexpansion = θ0[η] + ϵ θ1[η] + ϵ^2 θ2[η]+ ϵ^3 θ3[η]+ ϵ^4 θ4[η]

The nonlinear ODE is

ode =  ϵ θ[η]^2 -     D[θ[η],{η,2}]

Here is how to substitute for the dependent variable and any derivatives

result1 =ode/.{θ[η] ->perturbexpansion,Derivative[a1_][θ][η]:>  (∂_({η,a1}))perturbexpansion}

We need to collect the separate powers of ϵ.  We need each power of ϵ to be separately equal to 0 otherwise, the results will be correct, or not correct depending on a specific value of ϵ.  We need an answer for all ϵ

result2 = Collect[ExpandAll[result1],ϵ]

order0 = Coefficient[result2,ϵ,0]

order1 =Coefficient[result2,ϵ,1]

order2 = Coefficient[result2,ϵ,2]

order3 = Coefficient[result2,ϵ,3]

order4 = Coefficient[result2,ϵ,4]

Let’s get the anwers

order0answer =DSolve[{order0 ==0, θ0[0] == 1, θ0'[1] == 0}, θ0[η],η]

order1answer = DSolve[{((order1 )/. order0answer[[1]]) ==0, θ1[0] == 0, θ1'[1] == 0}, θ1[η],η]

order2answer = DSolve[{(((order2 )/. order1answer[[1]])/.order0answer[[1]]) ==0, θ2[0] == 0, θ2'[1] == 0}, θ2[η],η]

order3answer = DSolve[{((((order3 )/. order2answer[[1]])/.order1answer[[1]])/.order0answer[[1]]) ==0, θ3[0] == 0, θ3'[1] == 0}, θ3[η],η]

order4answer = DSolve[{(((((order4 )/. order3answer[[1]])/.order2answer[[1]])/.order1answer[[1]])/.order0answer[[1]]) ==0, θ4[0] == 0, θ4'[1] == 0}, θ4[η],η]

fourtermsanswer = ((((perturbexpansion/.(order4answer[[1]])/. order3answer[[1]])/.order2answer[[1]])/.order1answer[[1]])/.order0answer)

perturbploto3 =Plot[threetermsanswer/.{ϵ->{.01,.02,.05,.2,.5}},{η,0,1},PlotStyle->Dashing[{.01,.01}]]

perturbploto4 =Plot[fourtermsanswer/.{ϵ->{.01,.02,.05,.2,.5}},{η,0,1},PlotStyle->Dashing[{.02,.005}]]

Compare the 3 term solution to the numerical result

numericalplot

Show[perturbploto3,numericalplot]

Show[perturbploto4,numericalplot]

So there is some improvement, but, you have to stop sometime!!