A very simple model works nearby the normal area of operation for an LED. Just the following very approachable linear function:
$$V_D=V_F+I_D\cdot R_F$$
Note that my \$V_F\$ isn't the same as the one you wrote in your question. This is the minimal value which occurs when \$I_D=0\:\textrm{A}\$. (Which is never useful, so again this is just a model value.)
Also, \$V_D=V - I_D\cdot R_1\$, so:
$$\begin{align*}
V - I_D\cdot R_1&=V_F+I_D\cdot R_F\\\\
&\therefore\\\\
I_D&=\frac{V-V_F}{R_1+R_F}
\end{align*}$$
Let's say that \$V_F=1.6\:\textrm{V}\$ and \$R_F=25\:\Omega\$. Then you'd find that you get \$2.1\:\textrm{V}\$ at \$20\:\textrm{mA}\$. (Which given only one data point is perhaps the best we can do for now.)
You can now compute the value of \$R_1\$ as a function of the desired current:
$$R_1=\frac{V-V_F}{I_D}-R_F$$
And that's pretty easy to apply. Going with the above parameters, I get \$R_1=\frac{5\:\textrm{V}-1.6\:\textrm{V}}{20\:\textrm{mA}}-25\:\Omega=145\:\Omega\$.
But now you actually have a function allowing you to plug in smaller or larger current values that will use the local slope of the LED operation to get a better estimate for your resistor. (Assuming you got anywhere close with the parameter values.)
You only need two measurements with the actual LED using two nearby but different resistor values in order to figure out, pretty accurately, what the parameter values should be near the average measured current for that experiment. So it isn't difficult to come by, just from a bench test.
The actual diode voltage vs diode current curve isn't linear. But for a narrow range around the designed operation point it does have a slope that doesn't vary much, which is why the above equation can work reasonably well if you stay near the recommended operating point.
If you want to support a much wider range of behavior things get a little more complicated. Let's start by just showing the basic circuit you are talking about:
simulate this circuit – Schematic created using CircuitLab
The LED diode follows the Shockley equation model pretty well, over a very wide range of behaviors. To compute the intersection of the resistor load line and the LED diode curve, let's go through the crazy math involved now.
The current based upon the Shockley equation is:
$$I_{D}=I_{S}\cdot\left(e^\frac{V_{D}}{n\cdot V_T}-1\right)$$
In the above equation, \$I_S\$ is the saturation current of the diode or LED (which is itself a strong function of temperature), \$n\$ is the emission coefficient, and \$V_T\$ is the thermal voltage (about \$26\:\textrm{mV}\$ at room temperature.) The first two are model parameters and the last one is a physical characteristic that arrives from the large population statistics of colliding, interacting matter.
It is also the case that the diode voltage is what's left over after the resistor drops its voltage, so recalling what I wrote earlier:
$$V_D=V - I_D\cdot R_1$$
Putting these two together gives:
$$I_{D}=I_{S}\cdot\left(e^\frac{V - I_D\cdot R_1}{n\cdot V_T}-1\right)$$
Note that \$I_D\$ is on both sides of the equation.
Solving this requires the LambertW function, which solves \$v=u e^u\$ for \$u\$, given \$v\$. So all we have to do is get things into that form:
$$\begin{align*}
I_{D}&=I_{S}\cdot\left(e^\frac{V - I_D\cdot R_1}{n\cdot V_T}-1\right)\tag{1}\\\\
\left(I_{D}+I_{S}\right)e^\frac{I_D\cdot R_1}{n\cdot V_T}&=I_{S}\:e^\frac{V}{n\cdot V_T}\tag{2}\\\\
\frac{\left(I_D+I_S\right)\cdot R_1}{n\cdot V_T}\:e^\frac{I_D\cdot R_1}{n\cdot V_T}&=\frac{I_S\cdot R_1}{n\cdot V_T}\:e^\frac{V}{n\cdot V_T}\tag{3}\\\\
\frac{\left(I_D+I_S\right)\cdot R_1}{n\cdot V_T}\:e^\frac{\left(I_D+I_S\right)\cdot R_1}{n\cdot V_T}&=\frac{I_S\cdot R_1}{n\cdot V_T}\:e^\frac{V+I_S\cdot R_1}{n\cdot V_T}\tag{4}\\\\
&\therefore\\\\
\frac{\left(I_D+I_S\right)\cdot R_1}{n\cdot V_T}&=\operatorname{LambertW}\left(\frac{I_S\cdot R_1}{n\cdot V_T}\:e^\frac{V+I_S\cdot R_1}{n\cdot V_T}\right)\tag{5}\\\\
I_D &=\frac{n\cdot V_T}{R_1}\:\operatorname{LambertW}\left(\frac{I_S\cdot R_1}{n\cdot V_T}\:e^\frac{V+I_S\cdot R_1}{n\cdot V_T}\right)-I_S\tag{6}
\end{align*}$$
That's the actual mathematics involved. Usually, the value of \$I_S\$ is quite small, so the above can be simplified a bit:
$$\begin{align*}
I_D &\approx\frac{n\cdot V_T}{R_1}\:\operatorname{LambertW}\left(\frac{I_S\cdot R_1}{n\cdot V_T}\:e^\frac{V}{n\cdot V_T}\right)\tag{7}
\end{align*}$$
Of course, you need the model values for the diode. Different model values because this is based upon a different (and more complete) LED diode model.
Taking your example, I might guess that \$n=4\$ and \$I_S=35\:\textrm{pA}\$. Using the above equation I get \$I_D\approx 8.5\:\textrm{mA}\$. Note that this isn't the value you proposed. But I just made up some model values, too. And in reality, you'd need to hook things up just as you wrote and measure the results for either argument to be selected as closer. Who knows?
Absolutely no one I know ever does any of that, though. The above equation, if the model parameters are worked out and if the LED die temperature is maintained, will be very close to right over many, many orders of magnitude. It's surprisingly good over a large range of behaviors. But in practice for driving an LED, a designer doesn't need to go there.
There are lots of reasons why. The temperature of an LED is never really kept stable, in practice. And in any case, model parameters that work over a large range aren't needed because the LED is usually operated near its nominal current value. Besides, human perceptions of LED brightness are logarithmic and not particularly sensitive to modest differences in current (unless you are seeing two, side-by-side, I suppose.) So the whole point of the above exercise is more about being able to manipulate equations than being of any practical value for LEDs.
