1

The bvp system of ode's I am having for this problem is,

ode1 = f'''[y] + a * f'[y]*f'[y] - f[y]*f''[y] == 0
ode2 = T''[y] + b *f[y]*T'[y] == 0

with boundary conditions

bcs = {f[0] == 0, f'[0] == 1, f'[L] == 0, T[0] == 1, T[L] == 0};

where L=5. Taking some random values for a, b and c,NDSolve report stiffness. The expression which should take values from NDSolve is

SE = f'[y]^2 + b*T[y]^2 + c*f[y]*T[y]; (*at y = 0*)

I want to vary a$\in[0.1,0.5]$, b$\in[0.1,1]$ and c$\in[0,2]$ by 0.1 simultaneously for both the system and the secondary expression SE to create a data table for SE of the form,

a b c SE

Can someone please help?

Edit

The above system is just a test case. So ignore the convergence warnings. I am only interested in generating the table.

zhk
  • 11,939
  • 1
  • 22
  • 38
  • Can you add some background information for the PDEs? I believe it's not the first time I see this PDE set being asked here. Being aware of the background will also make the analysis easier. (By observing test = ParametricNDSolveValue[{ode1, bcs[[1 ;; 2]], f''[0] == ic} /. a -> 5/10, f, {y, 0, L}, ic]; Manipulate[ListLinePlot[test[ic]' // Quiet, PlotRange -> {{0, 5}, {-10, 10}}], {ic, -100, 100, 1/2}], I believe it's a weak solution problem. To overcome the stiffness, I guess a small modification is necessary for ode1, but finding it by pure math analysis is somewhat beyond my reach. ) – xzczd Dec 27 '16 at 07:28
  • @xzczd The idea behind these ODE's is from this post http://mathematica.stackexchange.com/questions/104170/heat-convection-differential-equations-from-1952-mathematica-fails-to-converg . My main interest is to get a table for the secondary equation SE. – zhk Dec 27 '16 at 07:37
  • @xzczd I have made a mistake in the first equation and now edit it. Sorry for the inconvenience. – zhk Dec 27 '16 at 07:39
  • Then when a=5/10, f'[L] seems to be always <0: te = ParametricNDSolveValue[{ode1, bcs[[;; 2]], f''[0] == ic} /. a -> 5/10, f, {y, 0, L}, ic]; Plot[te[ic]'[L], {ic, -3, 1}] – xzczd Dec 27 '16 at 08:31

1 Answers1

2

A straightforward usage of ParametricNDSolveValue and Table:

L = 5;
bcs = {f[0] == 0, f'[0] == 1, f'[L] == 0, T[0] == 1, T[L] == 0};
ode1 = f'''[y] + a*f'[y]*f'[y] - f[y]*f''[y] == 0;
ode2 = T''[y] + b f[y] T'[y] == 0;

pa = ParametricNDSolveValue[{ode1, ode2, bcs}, {f, T}, {y, 0, L}, {a, b}];

midtable = Flatten[
   Table[{a -> aval, b -> bval, {f, T} -> pa[aval, bval] // Thread} // Flatten, {aval, 
     1/10, 5/10, 1/10}, {bval, 1/10, 1, 1/10}], 1];

SE = f'[y]^2 + b*T[y]^2 + c*f[y]*T[y] /. y -> 0;
rst = Table[{a, b, c, SE} /. midtable // Evaluate, {c, 0, 2, 1/10}]
xzczd
  • 65,995
  • 9
  • 163
  • 468
  • The output of the Table is not elegant. It should be rows and columns not a list of data. – zhk Dec 30 '16 at 06:55
  • 3
    @MMM Flatten[rst, 1] // TableForm. Well, with all due respect, you should put a bit more effort into learning the core language of Mathematica. – xzczd Dec 30 '16 at 06:58