# NDSolve Plotting issue

I am trying to solve a system of ODEs with one extra boundary condition.

eq1 = f'''[y] + f[y] f''[y] - m f'[y] - f'[y]^2 - s (f'[y] + y/2 f''[y]) == 0;
eq2 = g''[y] - (s/2(3 g[y] + y g'[y]) + 2 g[y] f'[y] - g'[y] f[y]) + f''[y]^2 == 0;


With boundary conditions

bc = {f == 0, f' == 1, g == 1, f''[b] == 0, g'[b] == 0};


With the arbitrary trial values, the integrated solutions will generally not satisfy the outer boundary conditions.

f''[b] == 0 and f[b] == s b/2


I need to find a way to adjust the two trial values such that the numerical solution eventually matches the required boundary conditions and satisfy the relation f[b] = s b/2.

Here are my questions

• If I fix m (say = 1), then how I can plot b vs s?
• If I fix s (say = 1.2 and 0.8), then how I can plot f'[b] vs m?
• If I fix s (say = 1.2 and 0.8), then how I can plot f'' vs m?

What value b will take?

Here is my try, but it's not working

g[b_?NumericQ, m_?NumericQ, s_?NumericQ] :=
f[b] /.
NDSolve[{eq1, eq2, f == 0, f' == 1, g == 1, f''[b] == 0, g'[b] == 0},
{f, g}, {y, 0, b}]

s2[b_, s_] := FindRoot[g[b, s, m], {b, 0}]


Paper and Pen implementation of Shooting method As it is mentioned by @Peltio, here is the possible Shooting method.

Note $theta$ is g and $\eta$ is y. The output should look like this

1. b vs s 1. f'[b] vs m 1. f'' vs m • Please post your DSolve or NDSolve code ... and remember to start your own symbols with lowercase chars. – Dr. belisarius Nov 18 '13 at 1:12
• Something looks wrong: Either s2 should be a function of S and M, or the FindRoot variable should be M. – Michael E2 Nov 18 '13 at 1:38
• What is the relation of f to the other variables? Is it typo for F – m_goldberg Nov 18 '13 at 2:16
• Your equations aren't equations. When you have an equation, you need to use the equality operator, which is written == in Mathematica, rather than the assignment operator =. You should change both the ODEs and the boundary conditions. – Pillsy Nov 18 '13 at 2:52
• Thx I have removed the typo. – zhk Nov 18 '13 at 3:08

The following does work in principle but gives a lot of warnings and I don't think you should trust the numbers it gives, although they roughly seem to reproduce the first plot (which is the only I checked but don't know the m it was produced with). I have tried to keep this simple but robust and try to show some techniques that I found helpful in the past. So I hope it can serve as a guide about how to access similar problems.

# Avoid Redefinitions

to keep things simple I just do a

ClearAll[f, θ, η, f0, f1, f2, θ0, θ1, s, m]


in the beginning. If you plan to pack all that into a function I would certainly make all those local variables of a scoping construct, e.g. Module...

# Derive Equations

One of the probably underrated features of Mathematica is that you can input in almost exactly the form you find it in textbooks or publications. It will help a lot to avoid unnecessary typos if one makes good use of that. Thus I even sticked to the greek letters (which now even look nice here and can be copied to Mathematica thanks to shrx editing). To check with the original work you want to convert the results to StandardForm within Mathematica, e.g. via Cell->Convert To->StandardForm. It probably is worth noting that I detected such a type when trying this out before posting this, that typo is now fixed (missing (f1[η])^2 in first equation). I also let Mathematica do the simple algebra, which helps to avoid unneccessary typing and many of those annoying errors that people mentioned in the comments to your question. Here I was using this to define the equations:

eqns[m_, s_] = Function[
{f, θ, η},
Evaluate[{
f2'[η] == s*(f1[η] + η/2 f2[η]) ] + (f1[η])^2 - f0[η]*f2[η] + m*f1[η],
θ1'[η] == (s/2 (3*θ0[η] + η*θ1[η]) + 2*θ0[η] f1[η] - θ1[η]*f0[η]) - (f2[η])^2
} /. {
f0 -> f, f1 -> Derivative[f], f2 -> Derivative[f],
θ0 -> θ, θ1 -> Derivative[θ]
}
]
]

bc[β_] = Function[{f, θ, η},
Evaluate[{
f == 0, f' == 1, θ == 1,
f''[β] == 0, θ'[β] == 0
}]
]


note that by using Derivative one can write those replacments in a very elegant and robust way. I also consider it good practive to make the parameters function arguments in such cases.

# Find Solution

I split this in two parts, first is to solve the boundary condition problem. In newer versions of Mathematica NDSolve has a shooting method built in and will choose that automatically if given according boundary conditions. Thus this will give you a solution of f in the form of an InterpolatingFunction:

fSolution[β_?NumericQ, m_?NumericQ, s_?NumericQ] :=
First[NDSolveValue[
Evaluate[Flatten[{eqns[m, s][f, θ, η], bc[β][f, θ, η]}]],
{f}, {η, 0, β}
]]


This is a good stage where to test whether everything works as desired, e.g. by evaluating:

fSolution[1, 1, 1]


you would typically also investigate whether that solution fullfills the boundary conditions and in general looks like an appropriate solution...

Having this, we can now try to find that β which will fullfill the additional condition:

βSolution[m_?NumericQ, s_?NumericQ] := Module[{β},
β /. First[FindRoot[fSolution[β, m, s][β] == (s*β)/2, {β, 0.01, 10}]]
]


the initial values of FindRoot here are adjusted for the following table/plot, you probably might need to do something smarter to make this work for aribtrary cases...

# Calculate and Plot values

As the above is using two nested FindRoot (one with two variables implicitly by NDSolves shooting method) and an NDSolve within that calculating a solution isn't very fast. In such cases it often makes much more sense to decide at which points one wants to evaluate and calculate them instead of letting Plot do it's magic which in general will lead to a lot more evaluations. So here is quite rough approximation which is mainly good to check whether the above does what it should, for a nice plot you'd want to use some more points:

βOverS = Module[{x},
Table[x = {s, βSolution[1, s]}; Print[x]; x, {s, 0.2, 1, 0.1}]
]

ListLinePlot[βOverS, PlotRange -> {0, 10}, Frame -> True]


I have introduced a very crude monitoring of the progress. As mentioned above both FindRoot and NDSolve give some warnings about failing to find a good result and it will not even work for smaller values of s. I guess that you could get better results and get them faster with adjusting the various methods of the two functions. But that is a wide field and probably would justify a new question which should only concentrate on that aspect...

first, to make the following more comfortable one might want to memoize all beta-Values that were calculated. This can be done by redefining βSolution like this:

βSolution[m_?NumericQ, s_?NumericQ] := βSolution[m,s] = Module[{β},
Print[βSolution->{m,s}];
β /. First[FindRoot[fSolution[β, m, s][β] == (s*β)/2, {β, 0.01, 10}]]
]


I also added a Print so one can immediately see what is currently calculated. With these preparations one can then (almost) reproduce the other two plots like this:

fPrimeOverM[s_] := Module[{β},
Table[
β = βSolution[m, s];
{m, Derivative[fSolution[β, m, s]][β]},
{m, 0, 8, 0.5}
]
];

ListLinePlot[{fPrimeOverM[1.2], fPrimeOverM[0.8]},
Frame -> True,PlotStyle -> {Red, Black},
PlotRange -> {0, 0.5},AspectRatio -> 0.8
]


and:

fSecondOverM[s_] := Module[{β},
Table[
β = βSolution[m, s];
{m, -Derivative[fSolution[β, m, s]]},
{m, 0, 8, 0.5}
]
]

ListLinePlot[Evaluate[{fSecondOverM[1.2], fSecondOverM[0.8]}],
Frame -> True,
PlotStyle -> {Red, Black},
PlotRange -> {0, 3.5},
AspectRatio -> 0.8]


On my machine that gives larger discrepancies for small values of M, and that is also where FindRoot and NDSolve give a lot of warnings which indicates that these should be taken seriously. To really get reliable results you would definitely have to improve the settings of especially NDSolve (and of course recheck all definitions, including potential typos in the original work).

• Thx dear @Albert Retey, I have checked it and its working quit good for the first plot and some other similar plots like the first one. Any hint for the rest of the two plots? – zhk Nov 21 '13 at 3:16
• @MMM: using βSolution you can get the β for a given set of s,m. Inserting that in fSolution you can get the solution for f for those s,m in the form of an InterpolatingFunction. With that you can get the derivative and evaluate at β. Doing that for a range of m-Values and a fixed s will give you one of the lines of the second plot. Similary you can get the second derivative at 0, generate a list of values and plot that. I don't have time to do this right now but maybe later or on the weekend... – Albert Retey Nov 21 '13 at 10:36
• Thx alot @Albert Retey for such a nice and time consuming effort. Just one last thing, How we can extract the data from these plots to a .txt file? – zhk Nov 22 '13 at 14:39
• @MMM: Export["path/to/file.txt",βOverS,"Text"] – Albert Retey Nov 22 '13 at 15:09
• Dear @Albert Retey but the Export give the output with {,}? – zhk Nov 22 '13 at 15:42

This is not the answer to your problem, but perhaps this can help you find a solution to your plotting problem. I would start with defining the equations in terms of the parameters s and m:

eq1[m_, s_] = f'''[y] + f[y] f''[y] - m f'[y] - f'[y]^2 - s (f'[y] + y/2 f''[y]) == 0;
eq2[m_, s_] = g''[y] - (s/2(3 g[y] + y g'[y]) + 2 g[y] f'[y] - g'[y] f[y]) + f''[y]^2 == 0;


Then I will solve the systems of ODEs with the condition at one extreme of the interval.

vals[m_, s_, b_] := (
{f1, g1} = ({f, g} /. NDSolve[{eq1[m, s], eq2[m, s],
f == 0, f' == 1, f'' == -2,
g == 1, g' == -.2
}, {f, g}, {y, 0, b}][]);
{f1[b], g1[b]})
(*You can choose to return the derivatives if you wish*)


Mind you: this is not your system of ODEs+BCs. I had to make up further initial conditions for f and g. Your problem is more complicated in that you need to impose conditions at three points. (If m and s were known, this would be a straightforward problem of shooting to a fitting point - you would solve it using the extremes as initial/final condition and use the solution to compute the value at the same common point).

Now, if you specify one of the three values, you can plot the solutions in terms of the remaining twos. For example, suppose the middle point is b = 2, then the values of f and g in 2 can be seen as functions of m and s:

{Plot3D[vals[m, s, 2][], {m, -3, 3}, {s, -3, 3}],
Plot3D[vals[m, s, 2][], {m, -3, 3}, {s, -3, 3}]}


Incidentally, if you can get the solution to your system, this will allow you to do what you ask. Just specify the fixed value and plot vals[] - or just vals if you choose to return a single function's value - versus the other two. Of course vals can be made to return f, f' or f''.

• You misstyped the term g[y] f''[y] term in eq1 which should be f[y] f''[y]. – zhk Nov 20 '13 at 6:41