2
$\begingroup$

So I wish to solve the following $i$ equations for $\eta_i$

$$p(x)+\rho_i z_i + \frac{(u_id)^2}{2(\eta_{i-1}(x)-\eta_i(x))^2}+\rho_i \eta_i(x) = B_i$$ where we have the following definitions and $x$ varying between 0 and 1.

$p(x) = 1-x \\ B_i = p(0)+2z_i\rho_i + \frac{u_i^2}{2} \\ d = 0.1 \\ z = \{1,0.9,0.8,...,0.1,0 \} \\ \rho = \{0,0.1,0.2,...,0.9,1\} \\ u = \{1,1,1,..,...,1\} \\ \eta_1 = 1$

You will notice $z$ and $\rho$ are very rough approximations to linear lines and $d$ the discretisation is large. I have the following code to solve this system, which is already very slow. I feel I am not writing the most efficient code or making use of Mathematica properly. I wish to decrease the size of $d$ and have a much faster code. I think it is Solve which is slowing things down.

ρ = Table[0.1 n, {n, 0, 10}]
z = Table[1 - 0.1 n, {n, 0, 10}]
u = Table[1, {n, 0, 10}]
η = Table[y, {n, 0, 10}]
η[[1]] =  1
B = Table[1 + 2 z[[n]] ρ[[n]] + u[[n]]^2/2, {n, 1, 10}]

Do[η[[i]] = 
  Solve[1 - x + ρ[[i]]*z[[i]] + (d*u[[i]])^2/2(η[[i - 1]] - η[[
        i]])^2 + ρ[[i]]*η[[i]] == B[[i]], y], {i, 2, 10}]

Any pointers on how to use Mathematica properly to solve this system quickly?

EDIT

So I have modified the code so that it now solves quickly for each η[[i]], however the problem is now plotting the results, because there is a nested Root in each η[[i]] there is an error regarding the equation not being well formed? Anyway the new fast code for calculating the η[[i]]

y /. Normal@
  Simplify@Solve[1 - x + ρ[[i]]*z[[i]] + (d*u[[i]])^2/(
       2 (η[[i - 1]] - y)^2) + ρ[[i]]*y == B[[i]], y, Reals][[1]]

Then drop the output from this into

Do[η[[i]] = 
   Root[d^2 u[[i]]^2 + 2 η[[-1 + i]]^2 - 2 x η[[-1 + i]]^2 - 2 B[[i]] η[[-1 + i]]^2 + 
      2 z[[i]] η[[-1 + i]]^2 ρ[[i]] + (-4 η[[-1 + i]] + 4 x η[[-1 + i]] + 
         4 B[[i]] η[[-1 + i]] - 4 z[[i]] η[[-1 + i]] ρ[[i]] + 
         2 η[[-1 + i]]^2 ρ[[i]]) #1 + (2 - 2 x - 2 B[[i]] + 2 z[[i]] ρ[[i]] - 
         4 η[[-1 + i]] ρ[[i]]) #1^2 + 2 ρ[[i]] #1^3 &, 1], {i, 2, 10}];

For some reason it doesn't run if you the y/.Normal line within the Do loop. Doesn't work for a Table either. Anyway, if you then try and plot this it runs forever, but just plotting the first 3 gives what looks to be the correct result, however it does come with a "not well formed polynomial" error.

Plot[{η[[1]], η[[2]], η[[3]]}, {x, 0, 1}]

enter image description here

$\endgroup$
  • $\begingroup$ Why doesn't $x$ have an associated numerical value? $\endgroup$ – J. M. will be back soon Jan 18 '17 at 12:12
  • $\begingroup$ $x$ varies between 0 and 1 $\endgroup$ – M Horsley Jan 18 '17 at 12:28
  • 1
    $\begingroup$ Your code contains Eta[[i]], while Eta[[2]] is unknown. $\endgroup$ – Feyre Jan 18 '17 at 13:05
  • $\begingroup$ @Feyre I don't understand what you mean. I'm just setting \Eta[[i]] to the solution iteratively $\endgroup$ – M Horsley Jan 18 '17 at 13:56
  • $\begingroup$ @Feyre It's solving for \[Eta][[2]] ... If you notice \Eta[[2]] == y $\endgroup$ – M Horsley Jan 18 '17 at 14:33
2
$\begingroup$

Plotting η[[i]] is so slow, because it involves numerically evaluating Root functions nested i deep. However, the time can be reduced to almost nothing by using ListLinePlot instead of Plot and evaluating η[[i]] numerically within the final Do of the question.

a = Table[1, {i, 1, 10}, {j, 1, 11}];
Do[a[[i, j]] = η[[i]] = (Root[d^2 u[[i]]^2 + 2 η[[-1 + i]]^2 - 2 x η[[-1 + i]]^2 - 
    2 B[[i]] η[[-1 + i]]^2 + 2 z[[i]] η[[-1 + i]]^2 ρ[[i]] + (-4 η[[-1 + i]] + 
    4 x η[[-1 + i]] + 4 B[[i]] η[[-1 + i]] - 4 z[[i]] η[[-1 + i]] ρ[[i]] + 
    2 η[[-1 + i]]^2 ρ[[i]]) #1 + (2 - 2 x - 2 B[[i]] + 2 z[[i]] ρ[[i]] - 
    4 η[[-1 + i]] ρ[[i]]) #1^2 + 2 ρ[[i]] #1^3 &, 1]) /. x -> .1 (j - 1), 
    {j, 1, 11},{i, 2, 10}]
ListLinePlot[a, DataRange -> {0, 1}]

enter image description here

$\endgroup$
  • $\begingroup$ This works excellently, thank you. $\endgroup$ – M Horsley Jan 18 '17 at 23:32
  • $\begingroup$ @ojlm By the way, if you need a function that produces numbers quickly rather than a Do loop that produces a table of numbers, Fold should work well. $\endgroup$ – bbgodfrey Jan 18 '17 at 23:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.