英文:
Solving for a specific result to an equation with one variable
问题
如果我想解一个只有一个可变变量的方程,以获得我想要的结果,我该如何做呢?
最初我被建议使用二分法算法,但在进一步研究后,我发现它只适用于在Python中求解平方根方程。
例如:
rho1 = 0.177
Nh = 9.3
Nv = 128
Qh = 10
H = 1
beta_h1 = .0000000001
def beta_h_func(beta_h):
return beta_h1*min(Nv/(Qh*Nh), 1) + rho1*H*min(Nv/(Qh*Nh), 1)
while beta_h_func(beta_h1) != 1:
if beta_h_func(beta_h1) < 1:
beta_h1 += .000000001
beta_h_func(beta_h1)
if beta_h_func(beta_h1) > 1:
beta_h1 -= .000000001
beta_h_func(beta_h1)
其中 beta_h1
是可增加或减少的变量。但是,这样做只会导致一个无限的来回在1的上方和下方,永远不会结束。为了得到一个确定的结果1,我应该如何修改 while 循环呢?
英文:
If I want to solve for a specific result to an equation with one variable available to increase or decrease to solve for the result I want how could I go about doing it?
I was originally suggested the bisection algorithm but after researching that further I found it to only be applicable for square root equations in python.
For example
rho1 = 0.177
Nh = 9.3
Nv = 128
Qh = 10
H = 1
beta_h1 = .0000000001
def beta_h_func(beta_h):
return beta_h1*min(Nv/(Qh*Nh), 1) + rho1*H*min(Nv/(Qh*Nh), 1)
while beta_h_func(beta_h1) != 1:
if beta_h_func(beta_h1) < 1:
beta_h1 += .000000001
beta_h_func(beta_h1)
if beta_h_func(beta_h1) > 1:
beta_h1 -= .000000001
beta_h_func(beta_h1)
Where beta_h1
is the variable available for increasing or decreasing.
But doing this only results in an infinite back and forth going above and below 1 never ending. What could I change the while function to in order to receive a decisive 1 as my result.
答案1
得分: 1
这个问题是线性的,可以很容易地手动解决。解析解为 beta_h1 = (1 - rho1*H*min(Nv/(Qh*Nh), 1))/min(Nv/(Qh*Nh), 1)
。在这个例子中,该值为0.823。
假设这只是一个例子,你实际的问题更加复杂,你所面临的是一个寻根问题。你可以使用 scipy.optimize.root
来找到解。Scipy包含其他寻根方法,包括一个二分法,可以用于这个问题(我不确定你为什么说二分法只适用于平方根问题,不管那是什么意思)。
from scipy.optimize import root
rho1 = 0.177
Nh = 9.3
Nv = 128
Qh = 10
H = 1
c = min(Nv/(Qh*Nh), 1)
def beta_h_func(beta_h1):
return c*(beta_h1 + rho1*H) - 1
sol = root(beta_h_func, 1)
print(sol)
输出:
message: The solution converged.
success: True
status: 1
fun: [ 0.000e+00]
x: [ 8.230e-01]
nfev: 4
fjac: [[-1.000e+00]]
r: [-1.000e+00]
qtf: [ 3.861e-13]
根据这个结果,你的根是 sol.x[0]
(索引给出的是 float
而不是 numpy 数组)。
对你的代码有几点评论。
-
由于你使用浮点数并检查相等性,你的
while
循环可能会永远运行。在比较浮点数时,应该检查它们是否接近(使用一些容差),而不是完全相等。 -
在更新
beta_h1
后评估函数是无用的,因为这些结果没有被使用。 -
你的函数参数是
beta_h
而不是beta_h1
。 -
为了提高速度,你可以预先计算
min(...)
,因为其中的值在整个问题中不会改变(至少对于你分享的代码来说)。我在上面的代码中这样做了。
英文:
This problem is linear and can easily be solved by hand. The analytical solution is beta_h1 = (1 - rho1*H*min(Nv/(Qh*Nh), 1))/min(Nv/(Qh*Nh), 1)
. In this case, that value is 0.823.
Assuming this is just an example and your actual problem is more complicated, what you have is a root-finding problem. For that you can use scipy.optimize.root
to find the solution. Scipy contains other root-finding methods, including a bisection method that will work for this problem (I'm not sure why you said the bisection method only works for square root problems, whatever that means).
from scipy.optimize import root
rho1 = 0.177
Nh = 9.3
Nv = 128
Qh = 10
H = 1
c = min(Nv/(Qh*Nh), 1)
def beta_h_func(beta_h1):
return c*(beta_h1 + rho1*H) - 1
sol = root(beta_h_func, 1)
print(sol)
Output:
message: The solution converged.
success: True
status: 1
fun: [ 0.000e+00]
x: [ 8.230e-01]
nfev: 4
fjac: [[-1.000e+00]]
r: [-1.000e+00]
qtf: [ 3.861e-13]
Taking that result, your root is given as sol.x[0]
(indexing gives you the float
rather than the numpy array).
A couple of comments on your code.
- Your
while
loop will likely run forever since you are working with floating point numbers and are checking equality. When comparing floating point numbers, check if they're close (with some tolerance), not exactly equal. - The evaluation of the function after updating
beta_h1
is useless since nothing happens with those results. - The argument to your function is
beta_h
instead ofbeta_h1
. - For speed, you can precompute
min(...)
since the values inside do not change throughout your problem (at least not for the code you shared). I did this in my code above.
答案2
得分: 1
假设你的函数在替换正确答案后确实等于1,
>>> def beta_h_func(beta_h):
... return beta_h*min(Nv/(Qh*Nh), 1) + rho1*H*min(Nv/(Qh*Nh), 1)
...
>>> beta_h_func(.823)==1
True
你的循环逻辑是否正确?让我们在一个更简单的函数上进行测试:f(x) = x - 1
。
f = lambda x: x - 1
def go(step):
x = 0
dir = None
while (y:=f(x)) != 1:
if y < 1:
if dir is None:
dir = 1
elif dir == -1:
return 'looping', x, x+step
x += step
elif y > 1:
if dir is None:
dir = -1
elif dir == 1:
return 'looping', x-step, x
x -= step
return(x)
这个 go
函数接受一个步长,并返回使得 f(x)
等于1的 x
的值。它还具有一个安全功能,如果检测到我们在 x
的值之间循环(就像你遇到的情况),它会停止。让我们进行测试:
>>> go(1)
2
>>> go(.5)
2.0
>>> go(.25)
2.0
所以逻辑是正确的。让我们尝试一个更小的步长,比如0.1
>>> go(.1)
('looping', 1.9000000000000004, 2.0000000000000004)
尽管逻辑在精确算术中是正确的,但是当处理浮点数时,不能保证步长的总和将使你达到函数达到所需值的单个值。在我的例子中,从0开始以0.1为步长,我们到达了略大于1.9和略大于2的值,但从未到达2。你自己承认,你的函数也遭受了同样的命运。(它还需要近100亿步才能进入循环。)
值得花时间了解如何实现二分法,例如参考这里。
高级主题
如果你可以用整数输入来编写函数,甚至可以使用Python的内置二分法例程。例如,假设我们想找到方程 x**2 - 70
的解。让 x = i/prec
,然后定义一个函数类,它具有 __getitem__
方法,用于返回给定 i
的函数值,并允许你指定 prec
的值:
>>> class f(object):
... def __init__(self, prec): self.prec = prec
... def __getitem__(self, i):
... x = i/self.prec
... return x**2-70
...
>>> prec = 1000
>>> f(prec)[1*prec] # 1**2 - 70
-69.0
现在让 bisect
例程找到第一个使得函数值大于等于0的 x
的值:
>>> from bisect import bisect
>>> prec = 10000; bisect(f(prec), 0, 8*prec, 10*prec)/prec # 最接近1/1000
8.3667
>>> _**2
70.00166888999999
>>> prec = 2; bisect(f(prec), 0, 8*prec, 10*prec)/prec # 最接近一半
8.5
>>> _**2
72.25
>>> prec = 2**10; bisect(f(prec), 0, 8*prec, 10*prec)/prec # 最接近2**-10
8.3671875
>>> _**2
70.00982666015625
英文:
Assuming your function will actually equal 1 if you substitute in the correct answer,
>>> def beta_h_func(beta_h):
... return beta_h*min(Nv/(Qh*Nh), 1) + rho1*H*min(Nv/(Qh*Nh), 1)
...
>>> beta_h_func(.823)==1
True
is your loop logic sound? Let's test it on a simpler function: f(x) = x - 1
.
f = lambda x: x - 1
def go(step):
x = 0
dir = None
while (y:=f(x)) != 1:
if y < 1:
if dir is None:
dir = 1
elif dir == -1:
return 'looping', x, x+step
x += step
elif y > 1:
if dir is None:
dir = -1
elif dir == 1:
return 'looping', x-step, x
x -= step
return(x)
This go
function will accept a step and return the value of x
at which f(x)
is 1. It also has a safety feature to stop if it detects that we are looping between values of x
(like you are experiencing). Let's test it:
>>> go(1)
2
>>> go(.5)
2.0
>>> go(.25)
2.0
So the logic is sound. Let's try a smaller step...like 0.1
>>> go(.1)
('looping', 1.9000000000000004, 2.0000000000000004)
Although the logic is sound for exact artihmetic, when you are dealing with floating point numbers you are not assured that the sum of steps will bring you to the single value at which the function will attain the desired value. In my example, stepping by 0.1 from 0 brought us to just over 1.9 and just over 2, but never to 2. Your function (by your own admission) suffers from the same fate. (It also suffers by taking almost 10 billion steps to get to the point at which it loops.)
It is worth your time to see how to implement bisection, e.g. see here.
Advanced Topic
You can even use the built-in bisection routine of Python if you can write your function in terms of integer input. For example, say we wanted to find the solution of x**2 - 70
. Let x = i/prec
then define a function class that has a __getitem__
method to return the value of the function for a given i
and also allows you to specify the value of prec
:
>>> class f(object):
... def __init__(self, prec): self.prec = prec
... def __getitem__(self, i):
... x = i/self.prec
... return x**2-70
...
>>> prec = 1000
>>> f(prec)[1*prec] # 1**2 - 70
-69.0
Now let the bisect
routine find the first value of x that gives a value of f that is >= 0:
>>> from bisect import bisect
>>> prec = 10000; bisect(f(prec), 0, 8*prec, 10*prec)/prec # nearest 1/1000
8.3667
>>> _**2
70.00166888999999
>>> prec = 2; bisect(f(prec), 0, 8*prec, 10*prec)/prec # nearest half
8.5
>>> _**2
72.25
>>> prec = 2**10; bisect(f(prec), 0, 8*prec, 10*prec)/prec # nearest 2**-10
8.3671875
>>> _**2
70.00982666015625
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论