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$
8
  • $\begingroup$ Why doesn't $x$ have an associated numerical value? $\endgroup$ Jan 18, 2017 at 12:12
  • $\begingroup$ $x$ varies between 0 and 1 $\endgroup$
    – mch56
    Jan 18, 2017 at 12:28
  • 1
    $\begingroup$ Your code contains Eta[[i]], while Eta[[2]] is unknown. $\endgroup$
    – Feyre
    Jan 18, 2017 at 13:05
  • $\begingroup$ @Feyre I don't understand what you mean. I'm just setting \Eta[[i]] to the solution iteratively $\endgroup$
    – mch56
    Jan 18, 2017 at 13:56
  • $\begingroup$ @Feyre It's solving for \[Eta][[2]] ... If you notice \Eta[[2]] == y $\endgroup$
    – mch56
    Jan 18, 2017 at 14:33

1 Answer 1

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$
2
  • $\begingroup$ This works excellently, thank you. $\endgroup$
    – mch56
    Jan 18, 2017 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, 2017 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.