Physics 229& 100 Homework #8
Name: Nasser Abbasi  

1.  a) A  physics problem leads to the differential equation y[x]-x α y^′[x]+y^′′[x]=0 . Use DETools to find the first 10 terms of the power series expansion of the solutions about x=0.  Find separate expressions for solutions that contain even and odd powers of x.

Remove["Global`*"]

Use the MMHTools function powerss to find the power series solution around zero

<<"MMHTools`DETools`"

ode = y''[x] - x α y '[x] + y[x] ;

maxPower = 10 ;

seriesSolution = powerss[ode, y, x, maxPower, 0, False]

Now find the series with terms that are ODD powers of x

oddSeriesSolution = Total[Coefficient[seriesSolution, Table[x^n, {n, 1, maxPower, 2}]] * Table[x^n, {n, 1, maxPower, 2}]]

Lets look at the odd terms again from above,do not make any changes here,just look at coeff.or x powers to see the general term form

Cases[Collect[oddSeriesSolution, x], any_  x^b_.→ {x^b, any}] //TableForm

x a[1]
x^3 -a[1]/6 + 1/6 α a[1]
x^5 a[1]/120 - 1/30 α a[1] + 1/40 α^2 a[1]
x^7 -a[1]/5040 + 1/560 α a[1] - (23 α^2 a[1])/5040 + 1/336 α^3 a[1]
x^9 a[1]/362880 - (α a[1])/22680 + (43 α^2 a[1])/181440 - (11 α^3 a[1])/22680 + (α^4 a[1])/3456

Now find the series with terms that are even powers of x. We take x^0 as even power?

evenSeriesSolution = a[0] + evenSeriesSolution

Lets look at the even terms again from above,do not make any changes here,just look at coeff.or x powers to see the general term form

Cases[Collect[evenSeriesSolution, x], any_  x^b_.→ {x^b, any}] //TableForm

x^2 -a[0]/2
x^4 a[0]/24 - 1/12 α a[0]
x^6 -a[0]/720 + 1/120 α a[0] - 1/90 α^2 a[0]
x^8 a[0]/40320 - (α a[0])/3360 + (11 α^2 a[0])/10080 - 1/840 α^3 a[0]
x^10 -a[0]/3628800 + (α a[0])/181440 - (α^2 a[0])/25920 + (α^3 a[0])/9072 - (α^4 a[0])/9450

b) For certain values of α, the solutions are polynomials. Write a formula for the values of α for which polynomial solutions exist.

Since the ODE=0, then all the terms for the odd power of x must add to zero, and all the terms of even powers of x must add to zero. we do not want a[0] nor a[1] to be zero, since these are our free parameters, so it must be that the terms involving α sum to zero for each term.

For the even series, from above, I solve for α to see the pattern of α values which makes the coeff of each powers of x zero

Solve[1/24 - 1/12α == 0, α]

Solve[-1/720 + 1/120α - 1/90α^2 == 0, α]

Solve[1/40320 - α /3360 + (11 α^2 )/10080 - 1/840 α^3 == 0, α]

Solve[-1/3628800 + α /181440 - α^2 /25920 + α^3 /9072 - α^4 /9450 == 0, α]

{{α→1/2}}

{{α→1/4}, {α→1/2}}

{{α→1/6}, {α→1/4}, {α→1/2}}

{{α→1/8}, {α→1/6}, {α→1/4}, {α→1/2}}

So the same for the odd series, from above, I solve for α to see the pattern of α values which makes the coeff of each powers of x zero

Solve[-1/6 + 1/6 α == 0, α]

Solve[1/120 - 1/30 α + 1/40 α^2 == 0, α]

Solve[-1/5040 + 1/560 α - (23 α^2 )/5040 + 1/336 α^3 == 0, α]

Solve[1/362880 - α /22680 + (43 α^2 )/181440 - (11 α^3 )/22680 + α^4 /3456 == 0, α]

{{α→1}}

{{α→1/3}, {α→1}}

{{α→1/5}, {α→1/3}, {α→1}}

{{α→1/7}, {α→1/5}, {α→1/3}, {α→1}}

So,  even powers n of x, we need alpha to be 1/n,....,1/8,1/6,,1/4,1/2
So,  odd powers n of x, we need alpha to be 1/n,....,1/7,1/5,1/3,1

For the above values then a polynomial solution exist

c) Find the polynomials p1,p2,p3,p4 corresponding to the 4 largest values of α, and verify by explicit substitution that they satisfy the differential equation.

The 4 largest α values are shown above. Use these to find series solution for each specific α

p1 = (oddSeriesSolution/.α→ 1)

x a[1]

Now sub the above polynomial into the ODE. set α=1 also in the ODE itself

ode/.α→ 1/.y→Function[x, Evaluate[p1]]//Expand

0

p2 = (evenSeriesSolution/.α→ 1/2)

a[0] - 1/2 x^2 a[0]

Now sub the above polynomial into the ODE. set α=1/2 also in the ODE itself

ode/.α→ 1/2/.y→Function[x, Evaluate[p2]]//Expand

0

p3 = (oddSeriesSolution/.α→ 1/3)

x a[1] - 1/9 x^3 a[1]

Now sub the above polynomial into the ODE. set α=1/3 also in the ODE itself

ode/.α→ 1/3/.y→Function[x, Evaluate[p3]]//Expand

0

p4 = (evenSeriesSolution/.α→ 1/4)

a[0] - 1/2 x^2 a[0] + 1/48 x^4 a[0]

Now sub the above polynomial into the ODE. set α=1/4 also in the ODE itself

ode/.α→ 1/4/.y→Function[x, Evaluate[p4]]//Expand

0

d)Use DETools to  discuss the possible types of asymptotic behavior of the solutions to the differential equation as x→∞

Remove["Global`*"]

<<"MMHTools`DETools`"

First find what type of point we have at infinity

ode = y''[x] - x α y '[x] + y[x] ;

ClassifyODE[ode, y, x, ∞, False]

irregular singular point

Since the point x=∞ is an irreqular singular point, we can't use the method of replacing x->1/z and then do a power series at z=0 then switch back to x to do the leading terms analysis since power series expansion around an irreqular singular point will fail as illustrated below

newODE = ChangeVariable[ode, x→1/z, y, h, x, z]

h[z] + 2 z^3 h^′[z] + z α h^′[z] + z^4 h^′′[z]

Now the following command will FAIL, so I comment it out. Just for documenation is here.

(*frobenius[newODE, h, z, 0, 8, False] ; *)

The equation h[z] + 2 z^3 h^′[z] + z α h^′[z] + z^4 h^′′[z]   does not appear to have a regular singular point at z = 0

Now we do the standard analysis for finding asym behaviour for x=infinity when irrreqular singulr point as follows

trialX = ^(λ x^β) x^m ;

pits = ode/.y→Function[x, Evaluate[trialX]]//Simplify ;

poly = Expand[pits/trialX]

exps = uniqueexps[poly, x]

{-2, -2 + β, β, -2 + 2 β, 0}

Use MMHTools to do dominant balance analysis

dombalance[exps, β]

{2, 0}

So we have 2 consistant β value. We need to consider each in turn

poly/.β→2 ;

sp1 = Series[%, {x, ∞, 3}]

(-2 α λ + 4 λ^2)/(1/x)^2 + (1 - m α + 2 λ + 4 m λ) + (-m + m^2) (1/x)^2 + O[1/x]^4

Select[sp1[[3]], # =!= 0&]

{-2 α λ + 4 λ^2, 1 - m α + 2 λ + 4 m λ, -m + m^2}

set the first 2 to zero and solve

solλm = Solve[Thread[Take[%, 2] == 0], {λ, m}]

{{λ→0, m→1/α}, {λ→α/2, m→ (-1 - α)/α}}

Now obtain the 2 asym solutions

{as1, as2} = Thread[trialX/.solλm]/.β→2

{x^1/α, ^(x^2 α)/2 x^(-1 - α)/α}

Now we do the same for β=0

poly/.β→0 ;

sp1 = Series[%, {x, ∞, 5}]

(1 - m α) + (-m + m^2) (1/x)^2 + O[1/x]^6

Select[sp1[[3]], # =!= 0&]

{1 - m α, -m + m^2}

Since λ was eliminated when β=0 (see poly), then we can't solve for it and m as well. So I will not consider β=0 any more. Hence I stop here.

Conclusion: asymptotic behaviors of the soultion are combination of   x^1/α,  ^(x^2 α)/2 x^(-1 - α)/α
To continue further, I need to set a value for α (eigenvalue of the ODE) and get a general linear combination of the above 2 asym solution, as in  A x^1/α+B^(x^2 α)/2 x^(-1 - α)/αthen solve for A,B , then can plot the solution for large x.  I think this problem is only asking me to find the 2 asym. solutions, so I stop here. Otherwise I need to try different α's and for each solve for A,B and look at the solutions at large X. But asym. solutions should all behave the same way regadless of α value for large x.

Nice job

2. Consider the differential equation y'' [x]- y [x] x^(1/2)= 0 .
Determine the asymptotic behavior of its solutions.
a) Show that the point at infinity is an irregular singular point. Do this by explicitly performing the change of variable  z->1/x , and analyze the behavior at z=0.

Remove["Global`*"] ;

<<"MMHTools`DETools`"

ode = y''[x] - x^(1/2) y[x] ; (*ChangeVariable[ode, x→1/z, y, h, x, z] *)

newODE = ode/.y→Function[x, 1/x]/.x→z//Simplify

(2 - z^(5/2))/z^3

We see that for z = 0, we obtain 1/0, hence since z = 0 is the same as x = ∞, then we see that point at infinity is an irreqular singular point .

b) Use the MathematicaHandbook function  ClassifyODE to verify your answer.

ClassifyODE[ode, y, x, ∞, True]

 DE in standard form =  -(1/zv)^(9/2) yz[zv] + (2 yz^′[zv])/zv + yz^′′[zv]

      coeffs of derivatives of order<2 =  {-(1/zv)^(9/2), 2/zv}     analytic?    {False, False}

irregular singular point

c) Plug in a trial function of the form ^(λ x^β) x^m, and then divide thru by this expression to remove the exponential terms and produce a generalized polynomial in x. This represents the relative error of our asymptotic solution. We want to adjust the constants λ, β, and m to make this as small as possible when x is large.

Ok, now obtain trialx from trialZ and analyz what happens at z=0 (or now at x=∞)

trialX = ^(x^β λ) x^m

^(x^β λ) x^m

pits = ode/.y→Function[x, Evaluate[trialX]]//Simplify

^(x^β λ) x^(-2 + m) (m^2 - x^(5/2) + x^β (-1 + β) β λ + x^(2 β) β^2 λ^2 + m (-1 + 2 x^β β λ))

The above must vanish when z->0 (i.e. when x->∞)
The e^x terms can not vanish, so look for the other terms.

poly = pits/trialX//Expand

-m/x^2 + m^2/x^2 - x^(1/2) - x^(-2 + β) β λ + 2 m x^(-2 + β) β λ + x^(-2 + β) β^2 λ + x^(-2 + 2 β) β^2 λ^2

d) Make a list " by hand " of all the exponents which appear in the expression. The terms with the biggest exponent produce the biggest error in our asymptotic expansion. What is the highest constant power of x?  Use your freedom to pick β to produce another term with the same highest power of x. If there are several possibilities, make sure your choice is consistent.

RowBox[{RowBox[{We,  , want,  , to,  , adjust,  , Cell[m,λ,β], so,  , poly}], →, 0 at x→∞}]

Collect[poly, x^n_]

(-m + m^2)/x^2 - x^(1/2) + x^(-2 + 2 β) β^2 λ^2 + x^(-2 + β) (-β λ + 2 m β λ + β^2 λ)

Powers of x that appears are

                                          x^(-2),  x^1/2,  x^(-2 (1 - β)),  x^(-2 + β)

We want to choose β so that all the powers of x go to zero.  

the term x^(-2) we do not need to worry about since it goes to zero as x goes large. We just need to figure how to make Let look at how to make x^1/2 go to zero as x->∞

This then must be balanced by x to the powers of  -2(1-β)  and x to the power of   -2+β

solβ1 = Solve[-2 (1 - β) == 1/2, β]//First

{β→5/4}

solβ2 = Solve[-2 + β == 1/2, β]//First

{β→5/2}

Check each solution for consistancy with the second criterion.

{x^(-2 (1 - β)), x^(-2 + β)}/.solβ1

{x^(1/2), 1/x^(3/4)}

This is consistant since one power is positive and the second has negative power. Now we check the second solution

{x^(-2 (1 - β)), x^(-2 + β)}/.solβ2

{x^3, x^(1/2)}

Which is NOT consistant, since both has positive powers in x. Hence the solution is

                                                                          β→5/4

e) check your answer using the MathematicaHandbook function dombalance  

exps = uniqueexps[poly, x]

{-2, 1/2, -2 + β, -2 + 2 β}

dombalance[exps, β]

{5/4}

Verified ok

f)  Now use your freedom to choose λ so that  the coefficient of the highest power of x vanishes.

Now that we know what β should be, use it on poly

poly/.β→5/4

-m/x^2 + m^2/x^2 - x^(1/2) + (5 λ)/(16 x^(3/4)) + (5 m λ)/(2 x^(3/4)) + (25 x^(1/2) λ^2)/16

Collect[%, x]

(-m + m^2)/x^2 + ((5 λ)/16 + (5 m λ)/2)/x^(3/4) + x^(1/2) (-1 + (25 λ^2)/16)

Coefficient[%, Sqrt[x]]

-1 + (25 λ^2)/16

Solve[% == 0, λ]

{{λ→ -4/5}, {λ→4/5}}

g)  Now use your freedom to choose m so that the coefficient of the highest  remaining power of x vanishes, for each of the solutions to part f). Substitute the parameter values you have obtained into the original trial solution to get the leading order asymptotic behavior.

We had 2 solutions for λ in part (f) . Hence for λ→4/5 :

poly/.λ→4/5

-m/x^2 + m^2/x^2 - x^(1/2) - 4/5 x^(-2 + β) β + 8/5 m x^(-2 + β) β + 4/5 x^(-2 + β) β^2 + 16/25 x^(-2 + 2 β) β^2

%/.β→5/4

-m/x^2 + m^2/x^2 + 1/(4 x^(3/4)) + (2 m)/x^(3/4)

Collect[%, x]

(-m + m^2)/x^2 + (1/4 + 2 m)/x^(3/4)

Coefficient[%, 1/x^3/4]

1/4 + 2 m

Solve[% == 0, m]

{{m→ -1/8}}

Now for the second solution, for λ→ -4/5 :

poly/.λ→ -4/5

-m/x^2 + m^2/x^2 - x^(1/2) + 4/5 x^(-2 + β) β - 8/5 m x^(-2 + β) β - 4/5 x^(-2 + β) β^2 + 16/25 x^(-2 + 2 β) β^2

%/.β→5/4

-m/x^2 + m^2/x^2 - 1/(4 x^(3/4)) - (2 m)/x^(3/4)

Collect[%, x]

(-m + m^2)/x^2 + (-1/4 - 2 m)/x^(3/4)

Coefficient[%, 1/x^3/4]

-1/4 - 2 m

Solve[% == 0, m]

{{m→ -1/8}}

The same as before . Hence use  m→ -1/8, λ→4/5, λ→ -4/5, β→5/4

Now substitute these parameters into the trialX solution

as1 = trialX/.{m→ -1/8, λ→4/5, β→5/4}

^(4 x^(5/4))/5/x^(1/8)

as2 = trialX/.{m→ -1/8, λ→ -4/5, β→5/4}

^(-(4 x^(5/4))/5)/x^(1/8)

The above as1 and as2 are the leading order terms for the 2 asym. solutions

h) Check your answer using the MathematicaHandbook function asymptotic .

asymptotic[ode, y, x, 2, False]

{{^(-(4 x^(5/4))/5)/x^(1/8), 1 + 1881/(51200 x^(5/2)) - 9/(160 x^(5/4))}, {^(4 x^(5/4))/5/x^(1/8), 1 + 1881/(51200 x^(5/2)) + 9/(160 x^(5/4))}}

We see that the leading terms (which are the first term in each list) match result from part g).

i) Find the exact solutions to the differential equation using DSolve. Consider the special case  with y[0]==1,y'[0]=-1.Calculate the asymptotic behavior of this solution by forming a general linear combination of the asymptotic behaviors found in part g) and matching the value and the derivative of the  asymptotic and exact solutions at a point ( the exat value is not important, so pick x=5) Plot the eaxct solution and your asymptotic solution on the same graph with different colors.

sol = y[x]/.Flatten[DSolve[{ode == 0, y '[0] == 1, y[0] == 1}, y[x], x]]

1/5 (2^(3/5) 5^(2/5) x^(1/2) BesselI[2/5, (4 x^(5/4))/5] Gamma[2/5] + 2^(2/5) 5^(3/5) x^(1/2) BesselI[-2/5, (4 x^(5/4))/5] Gamma[3/5])

gen = A as1 + B as2

(B ^(-(4 x^(5/4))/5))/x^(1/8) + (A ^(4 x^(5/4))/5)/x^(1/8)

eq1 = (gen == sol)/.x→5 ;

eq2 = (D[gen, x] == D[sol, x])/.x→5 ;

{solA, solB} = {A, B}/.Flatten[Solve[{eq1, eq2}, {A, B}]]//N

{1.03914, 126.917}

So now that we have A and B,  we have the general asym. solution

gen = gen/.{A→solA, B→solB}

(126.917 ^(-(4 x^(5/4))/5))/x^(1/8) + (1.03914 ^(4 x^(5/4))/5)/x^(1/8)

Now compare the exact solution from DSolve with the asym. solution

Needs["Graphics`Colors`"]

[Graphics:HTMLFiles/index_206.gif]

We see that for large x, the asym. solution obtained is in very good agreement with Dsolve soltution The asym. solution for large x is in good agreement all the way down to x=4.

3. More about foxes and rabbits.  Assume that the fox and rabbit populations are governed by the differential equations  (d r[t])/(d t)= r[t] -r[t]^2/(2 f[t])  and (d f[t])/(d t)= 5 f[t] -5 f[t]^2 -(10 f[t] r[t])/(f[t] + .25) .  
a) The populations will be constant for conditions where r'[t] and f'[t] ==0. Find the stable poulations.

Remove["Global`*"]

rabbitODE = r '[t] == r[t] - r[t]^2/(2 f[t]) ;

foxODE = f '[t] == 5 f[t] - 5 f[t]^2 - (10 f[t] r[t])/(f[t] + 0.25) ;

To solve, use one of the above 2 rate equations. Set it to zero and solve for one population in terms of the second

f[t]/.Flatten[Solve[{r[t] - r[t]^2/(2 f[t]) == 0}, f[t]]]

r[t]/2

Hence, for stable conditions, the fox population will be 1/2 that of the rabbits . Now I will replace r[t] by 2 * f[t] in the Fox ODE and solve for the f[t] numerical valuae

foxODE/.r[t] →2 f[t]

f^′[t] == 5 f[t] - 5 f[t]^2 - (20 f[t]^2)/(0.25 + f[t])

%/.lhs_ == rhs_→rhs

5 f[t] - 5 f[t]^2 - (20 f[t]^2)/(0.25 + f[t])

Flatten[Solve[% == 0, f[t]]]

{f[t] → -3.32518, f[t] →0., f[t] →0.0751838}

From the above, since I am looking for fox population that is ≥ 0 to make any sense, then there is only only choice which is 0.075

foxPopulation = f[t]/.%[[3]]

0.0751838

and since rabiits are 2 * foxes then we have

rabbitsPopulation = 2 * %

0.150368

b) Show that there is a basin of attraction by solving the differential equations for an initial state slightly different from the equilibrium state, for example r[0]=.5, f[0]=0.3. Make  plots of the solution which show r[t] vs. t, f[t] vs t and a parametric plot in the r-t plane. Show that the population tends toward the equilibrium you found in part a.

sol = Flatten[NDSolve[{rabbitODE, foxODE, r[0] == 0.5, f[0] == 0.3}, {r[t], f[t]}, {t, 0, 20}]]

{r[t] →InterpolatingFunction[{{0., 20.}}, <>][t], f[t] →InterpolatingFunction[{{0., 20.}}, <>][t]}

ParametricPlot[{r[t], f[t]}/.sol, {t, 0, 20}, PlotRange→ {{0, 0.6}, {0, 0.4}}, AxesOrigin→ {0, 0}, AxesLabel→ {"rabbits\npopulation", "fox\npopulation"}] ;

[Graphics:HTMLFiles/index_231.gif]

We see from the plot above that equilibrium is at where fox population is half that of the rabbits, as found in part (a). To make it more clear, zoom in

sol = Flatten[NDSolve[{rabbitODE, foxODE, r[0] == 0.5, f[0] == 0.3}, {r[t], f[t]}, {t, 0, 40}]]

ParametricPlot[{r[t], f[t]}/.sol, {t, 0, 40}, PlotRange→ {{0.12, 0.16}, {0.05, 0.08}}, AxesLabel→ {"rabbits\npopulation", "fox\npopulation"}] ;

{r[t] →InterpolatingFunction[{{0., 40.}}, <>][t], f[t] →InterpolatingFunction[{{0., 40.}}, <>][t]}

[Graphics:HTMLFiles/index_235.gif]

4.a) At t=0 a  duck starts at (1,0) and paddles at unit velocity counterclockwise around the unit circle. At t=0, a dog initially at the origin swims with constant unit speed with velocity directed at the instantaneous position of the duck. Find the position of the dog as a function of time; the equations are non-linear, so  use NDSolve. Make a plot of the trajectory of the dog using ParametricPlot.  

Hint:

Remove["Global`*"] ;

duck = {Cos[t], Sin[t]} ;

vduck = D[duck, t] ;

dog = {xd[t], yd[t]} ;

dist = Sqrt[(duck - dog) . (duck - dog)] ;

vdog = (duck - dog)/dist ;

eqmo = Thread[D[dog, t] == vdog]

dogpath = NDSolve[Join[eqmo, {xd[0] == 0, yd[0] == 0}], {xd, yd}, {t, 0, 139}]

ParametricPlot[Evaluate[{xd[t], yd[t]}/.dogpath], {t, 0, 60}, AspectRatio→Automatic] ;

{xd^′[t] == (Cos[t] - xd[t])/((Cos[t] - xd[t])^2 + (Sin[t] - yd[t])^2)^(1/2), yd^′[t] == (Sin[t] - yd[t])/((Cos[t] - xd[t])^2 + (Sin[t] - yd[t])^2)^(1/2)}

{{xd→InterpolatingFunction[{{0., 139.}}, <>], yd→InterpolatingFunction[{{0., 139.}}, <>]}}

[Graphics:HTMLFiles/index_247.gif]

b) Does the dog catch the duck? Make a Log-Log plot of the distance between the dog and the duck to help decide this issue.

First make a plot of distance vs. time in linear scale to see what is going on

[Graphics:HTMLFiles/index_249.gif]

Now make a log-log plot. Similar to the above, but use Log to base 10 for the time and for the distance

<<Graphics`Graphics`

[Graphics:HTMLFiles/index_252.gif]

We see from the log-log, that distance never gets to zero. The dots get more and more dense as time passes (becuase dog keep changing its position more and more as it gets close to the duck).

c) Also make an animation which shows the motion of both the dog and the duck, using different symbols or colors to represent them. See   Plotting&GraphicsExamples.nb for a  guide on how to do the animation.

In this animation I scale everything by 10 to make it easier to see.

Needs["Graphics`Animation`"] (*first make a cicle plot to add it to animation*)

circleShape = ParametricPlot[10 {Cos[t], Sin[t]}, {t, 0, 2π}, AspectRatio→Automatic] ;

[Graphics:HTMLFiles/index_255.gif]

Needs["Graphics`Colors`"]

dogCoordinates = Flatten[Table[Evaluate[10 {xd[t], yd[t]}]/.dogpath, {t, 0, 7, .05}], 1] ;

dogObjects = Thread[ Circle[dogCoordinates, .3]] ;

duckCoordinates = Table[10 {Cos[t], Sin[t]}, {t, 0, 7, .05}] ;

duckObjects = Thread[ Circle[duckCoordinates, .3]] ;

graphicsTable = Table[{circleShape, Graphics[{{Red, dogObjects[[i]]} , duckObjects[[i]]} , PlotRange-> {{-15, 15}, {-15, 15}}   ]}, {i, 1, Length[dogObjects]}] ;

ShowAnimation[%] ;

SelectionMove[EvaluationNotebook[], All, GeneratedCell] ;

    FrontEndTokenExecute["CellGroup"] ;

      FrontEndTokenExecute["OpenCloseGroup"]

[Graphics:HTMLFiles/index_407.gif]

5. A non-linear boundary value problem. Make a plot of the solution of
y''[x] +x y'[x] +Sqrt[Abs[y[x]]]==0, with y[0]==0 and y[4]==0. To do this you will need to find a value of y'[0] that results in y[4]=0.
a) Make a function y4[yp_] that takes a value of y'[0] as input and returns the value of y[4]. Make a plot of this function over some reasonable range.

Remove["Global`*"] ;

ode = y''[x] + x y '[x] + Abs[y[x]]^(1/2) ;

y4[yp_ ? NumericQ] := Module[{sol}, y[4]/.NDSolve[{ode == 0, y '[0] == yp, y[0] == 0}, y, {x, 0, 10}]//First] ;

Plot[y4[initialSpeed], {initialSpeed, 0, 5}, AxesLabel→ {"initial speed", "position at x=4"}] ;

[Graphics:HTMLFiles/index_412.gif]

b) Use FindRoot to locate the value yp that does the job. Note that you will have to use FindRoot with 2 initial guesses (see Help -> FindRoot). Then plot the actual solution to the boundary value problem y[x].

requiredInitialSpeed = FindRoot[y4[v], {v, 0.5, 1.5}]

{v→0.787515}

sol = NDSolve[{ode == 0, y '[0] == v/.requiredInitialSpeed, y[0] == 0}, y, {x, 0, 15}]//First

{y→InterpolatingFunction[{{0., 15.}}, <>]}

Plot[y[x]/.sol, {x, 0, 15}, AxesLabel→ {"x", "y[x]"}, PlotLabel->"Solution using shooting method"] ;

[Graphics:HTMLFiles/index_418.gif]


Created by Mathematica  (January 30, 2007) Valid XHTML 1.1!