# Improving my method for solving iterative algebraic equations

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}]


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

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}]


• This works excellently, thank you. Commented Jan 18, 2017 at 23:32
• @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. Commented Jan 18, 2017 at 23:51