欧拉光束，求解python中的微分方程
Euler beam, solving differential equation in python
我必须求解欧拉伯努利微分光束方程：
w’’’’(x) = q(x)
和边界条件：
w(0) = w(l) = 0
和
w′′(0) = w′′(l) = 0
光束如下图所示：
beam
连续力q是2N/mm。
我必须使用拍摄方法和scipy.integrate.odeint()功能
我什至无法开始，因为我不明白如何将微分方程写成方程组
懂python边界条件微分方程求解的大神能帮忙吗！
谢谢:)
拍摄方法
要使用 shooting method 解决 fourth order ODE BVP 和 scipy.integrate.odeint() 你需要：
1.) 通过替换将 4th order ODE 分隔成 4 first order ODEs：
u = w
u1 = u' = w' # 1
u2 = u1' = w'' # 2
u3 = u2' = w''' # 3
u4 = u3' = w'''' = q # 4
2.) 创建一个 function 来执行 derivation logic 并将该函数连接到 integrate.odeint()，如下所示：
function calc(u, x , q)
{
return [u[1], u[2], u[3] , q]
}
w = integrate.odeint(calc, [w(0), guess, w''(0), guess], xList, args=(q,))
解释：
我们将 odeint() 的边界值条件发送给 x=0 ([w(0), w'(0) ,w''(0), w'''(0)])，它调用函数 calc，returns 导数是添加到 w 的当前状态。请注意，我们在输入已知的 w(0)=0 和 w''(0)=0.
时猜测 w'(0) 和 w'''(0) 的初始边界条件
将导数添加到 w 的当前状态如下：
# the current w(x) value is the previous value plus the current change of w in dx.
w(x) = w(x-dx) + dw/dx
# others are calculated the same
dw(x)/dx = dw(x-dx)/dx + d^2w(x)/dx^2
# etc.
这就是我们从 calc function 返回值 [u[1], u[2], u[3] , q] 而不是 [u[0], u[1], u[2] , u[3]] 的原因，因为 u[1] 是 first derivative 所以我们将它添加到w等
3.) 现在我们可以设置 shooting method。我们将 w'(0) 和 w'''(0) 的 different initial boundary values 发送到 odeint()，然后检查 returned w(x) profile 的最终结果以确定 w(L) 和w''(L) 达到 0（known 边界条件）。
拍摄方法程序：
# a function to return the derivatives of w
def returnDerivatives(u, x, q):
return [u[1], u[2], u[3], q]
# a shooting funtion which takes in two variables and returns a w(x) profile for x=[0,L]
def shoot(u2, u4):
# the number of x points to calculate integration -> determines the size of dx
# bigger number means more x's -> better precision -> longer execution time
xSteps = 1001
# length of the beam
L= 1.0 # 1m
xSpace = np.linspace(0, L, xSteps)
q = 0.02 # constant [N/m]
# integrate and return the profile of w(x) and it's derivatives, from x=0 to x=L
return odeint(returnDerivatives, [ 0, u2, 0, u4] , xSpace, args=(q,))
# the tolerance for our results.
tolerance = 0.01
# how many numbers to consider for u2 and u4 (the guess boundary conditions)
u2_u4_maxNumbers = 1327 # bigger number, better precision, slower program
# you can also divide into separate variables like u2_maxNum and u4_maxNum
# these are already tested numbers (the best results are somewhere in here)
u2Numbers = np.linspace(-0.1, 0.1, u2_u4_maxNumbers)
# the same as above
u4Numbers = np.linspace(-0.5, 0.5, u2_u4_maxNumbers)
# result list for extracted values of each w(x) profile => [u2Best, u4Best, w(L), w''(L)]
# which will help us determine if the w(x) profile is inside tolerance
resultList = []
# result list for each U (or w(x) profile) => [w(x), w'(x), w''(x), w'''(x)]
resultW = []
# start generating numbers for u2 and u4 and send them to odeint()
for u2 in u2Numbers:
for u4 in u4Numbers:
U = []
U = shoot(u2,u4)
# get only the last row of the profile to determine if it passes tolerance check
result = U[len(U)-1]
# only check w(L) == 0 and w''(L) == 0, as those are the known boundary cond.
if (abs(result[0]) < tolerance) and (abs(result[2]) < tolerance):
# if the result passed the tolerance check, extract some values from the
# last row of the w(x) profile which we will need later for comaprisons
resultList.append([u2, u4, result[0], result[2]])
# add the w(x) profile to the list of profiles that passed the tolerance
# Note: the order of resultList is the same as the order of resultW
resultW.append(U)
# go through the resultList (list of extracted values from last row of each w(x) profile)
for i in range(len(resultList)):
x = resultList[i]
# both boundary conditions are 0 for both w(L) and w''(L) so we will simply add
# the two absolute values to determine how much the sum differs from 0
y = abs(x[2]) + abs(x[3])
# if we've just started set the least difference to the current
if i == 0:
minNum = y # remember the smallest difference to 0
index = 0 # remember index of best profile
elif y < minNum:
# current sum of absolute values is smaller
minNum = y
index = i
# print out the integral for w(x) over the beam
sum = 0
for i in resultW[index]:
sum = sum + i[0]
print("The integral of w(x) over the beam is:")
print(sum/1001) # sum/xSteps
这输出：
The integral of w(x) over the beam is:
0.000135085272117
要打印出我们找到的 w(x) 的最佳配置文件：
print(resultW[index])
输出如下：
# w(x) w'(x) w''(x) w'''(x)
[[ 0.00000000e+00 7.54147813e-04 0.00000000e+00 -9.80392157e-03]
[ 7.54144825e-07 7.54142917e-04 -9.79392157e-06 -9.78392157e-03]
[ 1.50828005e-06 7.54128237e-04 -1.95678431e-05 -9.76392157e-03]
...,
[ -4.48774290e-05 -8.14851572e-04 1.75726275e-04 1.01560784e-02]
[ -4.56921910e-05 -8.14670764e-04 1.85892353e-04 1.01760784e-02]
[ -4.65067671e-05 -8.14479780e-04 1.96078431e-04 1.01960784e-02]]
为了仔细检查上面的结果，我们还将使用 numerical method. 求解 ODE
数值法
要使用numerical method解决问题我们首先需要solve the differential equations。我们将得到 four constants，我们需要在 boundary conditions 的帮助下找到它。 boundary conditions 将用于形成 system of equations 以帮助找到必要的 constants。
例如：
w’’’’(x) = q(x);
意味着我们有这个：
d^4(w(x))/dx^4 = q(x)
因为 q(x) 在积分后是常数，我们有：
d^3(w(x))/dx^3 = q(x)*x + C
再次整合后：
d^2(w(x))/dx^2 = q(x)*0.5*x^2 + C*x + D
再整合后：
dw(x)/dx = q(x)/6*x^3 + C*0.5*x^2 + D*x + E
最后的积分结果为：
w(x) = q(x)/24*x^4 + C/6*x^3 + D*0.5*x^2 + E*x + F
然后我们看一下 boundary conditions（现在我们有上面的 w''(x) 和 w(x) 的表达式）我们用它来制作 system of equations 来解决constants.
w''(0) => 0 = q(x)*0.5*0^2 + C*0 + D
w''(L) => 0 = q(x)*0.5*L^2 + C*L + D
这给了我们 constants:
D = 0 # from the first equation
C = - 0.01 * L # from the second (after inserting D=0)
对 w(0)=0 和 w(L)=0 重复相同的操作后，我们得到：
F = 0 # from first
E = 0.01/12.0 * L^3 # from second
现在，在我们有 solved the equation 并找到所有 integration constants 之后，我们可以为 numerical method.
编写程序
数值方法程序
我们将制作一个 FOR 循环，一次遍历每个 dx 的整个波束，并总结（积分）w(x)。
L = 1.0 # in meters
step = 1001.0 # how many steps to take (dx)
q = 0.02 # constant [N/m]
integralOfW = 0.0; # instead of w(0) enter the boundary condition value for w(0)
result = []
for i in range(int(L*step)):
x= i/step
w = (q/24.0*pow(x,4) - 0.02/12.0*pow(x,3) + 0.01/12*pow(L,3)*x)/step # current w fragment
# add up fragments of w for integral calculation
integralOfW += w
# add current value of w(x) to result list for plotting
result.append(w*step);
print("The integral of w(x) over the beam is:")
print(integralOfW)
输出：
The integral of w(x) over the beam is:
0.00016666652805511192
现在比较两种方法
射击法与数值法结果对比
w(x)在横梁上的积分：
Shooting method -> 0.000135085272117
Numerical method -> 0.00016666652805511192
这是一个很好的匹配，现在让我们看看图表：
从图中可以更明显地看出我们匹配得很好并且 shooting method 的结果是正确的。
要获得更好的 shooting method 结果，请将 xSteps 和 u2_u4_maxNumbers 增加到更大的数字，您还可以缩小 u2Numbers 和 u4Numbers到 same set size but a smaller interval（大约是上一个程序 运行 的最佳结果）。请记住，将 xSteps 和 u2_u4_maxNumbers 设置得太高会导致您的程序长时间 运行。
您需要将 ODE 转换为一阶系统，设置 u0=w 一个可能且通常使用的系统是
u0'=u1,
u1'=u2,
u2'=u3,
u3'=q(x)
这可以实现为
def ODEfunc(u,x): return [ u[1], u[2], u[3], q(x) ]
然后制作一个以实验初始条件和returns第二个边界条件的分量为准的函数
def shoot(u01, u03): return odeint(ODEfunc, [0, u01, 0, u03], [0, l])[-1,[0,2]]
现在你有一个包含两个分量的两个变量的函数，你需要用通常的方法求解这个 2x2 系统。由于系统是线性的，所以射击函数也是线性的，你只需要找到系数并求解得到的线性系统。
我必须求解欧拉伯努利微分光束方程：
w’’’’(x) = q(x)
和边界条件：
w(0) = w(l) = 0
和
w′′(0) = w′′(l) = 0
光束如下图所示：
beam
连续力q是2N/mm。
我必须使用拍摄方法和scipy.integrate.odeint()功能
我什至无法开始，因为我不明白如何将微分方程写成方程组
懂python边界条件微分方程求解的大神能帮忙吗！
谢谢:)
拍摄方法
要使用 shooting method 解决 fourth order ODE BVP 和 scipy.integrate.odeint() 你需要：
1.) 通过替换将 4th order ODE 分隔成 4 first order ODEs：
u = w
u1 = u' = w' # 1
u2 = u1' = w'' # 2
u3 = u2' = w''' # 3
u4 = u3' = w'''' = q # 4
2.) 创建一个 function 来执行 derivation logic 并将该函数连接到 integrate.odeint()，如下所示：
function calc(u, x , q)
{
return [u[1], u[2], u[3] , q]
}
w = integrate.odeint(calc, [w(0), guess, w''(0), guess], xList, args=(q,))
解释：
我们将 odeint() 的边界值条件发送给 x=0 ([w(0), w'(0) ,w''(0), w'''(0)])，它调用函数 calc，returns 导数是添加到 w 的当前状态。请注意，我们在输入已知的 w(0)=0 和 w''(0)=0.
w'(0) 和 w'''(0) 的初始边界条件
将导数添加到 w 的当前状态如下：
# the current w(x) value is the previous value plus the current change of w in dx.
w(x) = w(x-dx) + dw/dx
# others are calculated the same
dw(x)/dx = dw(x-dx)/dx + d^2w(x)/dx^2
# etc.
这就是我们从 calc function 返回值 [u[1], u[2], u[3] , q] 而不是 [u[0], u[1], u[2] , u[3]] 的原因，因为 u[1] 是 first derivative 所以我们将它添加到w等
3.) 现在我们可以设置 shooting method。我们将 w'(0) 和 w'''(0) 的 different initial boundary values 发送到 odeint()，然后检查 returned w(x) profile 的最终结果以确定 w(L) 和w''(L) 达到 0（known 边界条件）。
拍摄方法程序：
# a function to return the derivatives of w
def returnDerivatives(u, x, q):
return [u[1], u[2], u[3], q]
# a shooting funtion which takes in two variables and returns a w(x) profile for x=[0,L]
def shoot(u2, u4):
# the number of x points to calculate integration -> determines the size of dx
# bigger number means more x's -> better precision -> longer execution time
xSteps = 1001
# length of the beam
L= 1.0 # 1m
xSpace = np.linspace(0, L, xSteps)
q = 0.02 # constant [N/m]
# integrate and return the profile of w(x) and it's derivatives, from x=0 to x=L
return odeint(returnDerivatives, [ 0, u2, 0, u4] , xSpace, args=(q,))
# the tolerance for our results.
tolerance = 0.01
# how many numbers to consider for u2 and u4 (the guess boundary conditions)
u2_u4_maxNumbers = 1327 # bigger number, better precision, slower program
# you can also divide into separate variables like u2_maxNum and u4_maxNum
# these are already tested numbers (the best results are somewhere in here)
u2Numbers = np.linspace(-0.1, 0.1, u2_u4_maxNumbers)
# the same as above
u4Numbers = np.linspace(-0.5, 0.5, u2_u4_maxNumbers)
# result list for extracted values of each w(x) profile => [u2Best, u4Best, w(L), w''(L)]
# which will help us determine if the w(x) profile is inside tolerance
resultList = []
# result list for each U (or w(x) profile) => [w(x), w'(x), w''(x), w'''(x)]
resultW = []
# start generating numbers for u2 and u4 and send them to odeint()
for u2 in u2Numbers:
for u4 in u4Numbers:
U = []
U = shoot(u2,u4)
# get only the last row of the profile to determine if it passes tolerance check
result = U[len(U)-1]
# only check w(L) == 0 and w''(L) == 0, as those are the known boundary cond.
if (abs(result[0]) < tolerance) and (abs(result[2]) < tolerance):
# if the result passed the tolerance check, extract some values from the
# last row of the w(x) profile which we will need later for comaprisons
resultList.append([u2, u4, result[0], result[2]])
# add the w(x) profile to the list of profiles that passed the tolerance
# Note: the order of resultList is the same as the order of resultW
resultW.append(U)
# go through the resultList (list of extracted values from last row of each w(x) profile)
for i in range(len(resultList)):
x = resultList[i]
# both boundary conditions are 0 for both w(L) and w''(L) so we will simply add
# the two absolute values to determine how much the sum differs from 0
y = abs(x[2]) + abs(x[3])
# if we've just started set the least difference to the current
if i == 0:
minNum = y # remember the smallest difference to 0
index = 0 # remember index of best profile
elif y < minNum:
# current sum of absolute values is smaller
minNum = y
index = i
# print out the integral for w(x) over the beam
sum = 0
for i in resultW[index]:
sum = sum + i[0]
print("The integral of w(x) over the beam is:")
print(sum/1001) # sum/xSteps
这输出：
The integral of w(x) over the beam is:
0.000135085272117
要打印出我们找到的 w(x) 的最佳配置文件：
print(resultW[index])
输出如下：
# w(x) w'(x) w''(x) w'''(x)
[[ 0.00000000e+00 7.54147813e-04 0.00000000e+00 -9.80392157e-03]
[ 7.54144825e-07 7.54142917e-04 -9.79392157e-06 -9.78392157e-03]
[ 1.50828005e-06 7.54128237e-04 -1.95678431e-05 -9.76392157e-03]
...,
[ -4.48774290e-05 -8.14851572e-04 1.75726275e-04 1.01560784e-02]
[ -4.56921910e-05 -8.14670764e-04 1.85892353e-04 1.01760784e-02]
[ -4.65067671e-05 -8.14479780e-04 1.96078431e-04 1.01960784e-02]]
为了仔细检查上面的结果，我们还将使用
numerical method. 求解 ODE
数值法
要使用numerical method解决问题我们首先需要solve the differential equations。我们将得到 four constants，我们需要在 boundary conditions 的帮助下找到它。 boundary conditions 将用于形成 system of equations 以帮助找到必要的 constants。
例如：
w’’’’(x) = q(x);
意味着我们有这个：
d^4(w(x))/dx^4 = q(x)
因为 q(x) 在积分后是常数，我们有：
d^3(w(x))/dx^3 = q(x)*x + C
再次整合后：
d^2(w(x))/dx^2 = q(x)*0.5*x^2 + C*x + D
再整合后：
dw(x)/dx = q(x)/6*x^3 + C*0.5*x^2 + D*x + E
最后的积分结果为：
w(x) = q(x)/24*x^4 + C/6*x^3 + D*0.5*x^2 + E*x + F
然后我们看一下 boundary conditions（现在我们有上面的 w''(x) 和 w(x) 的表达式）我们用它来制作 system of equations 来解决constants.
w''(0) => 0 = q(x)*0.5*0^2 + C*0 + D
w''(L) => 0 = q(x)*0.5*L^2 + C*L + D
这给了我们 constants:
D = 0 # from the first equation
C = - 0.01 * L # from the second (after inserting D=0)
对 w(0)=0 和 w(L)=0 重复相同的操作后，我们得到：
F = 0 # from first
E = 0.01/12.0 * L^3 # from second
现在，在我们有 solved the equation 并找到所有 integration constants 之后，我们可以为 numerical method.
数值方法程序
我们将制作一个 FOR 循环，一次遍历每个 dx 的整个波束，并总结（积分）w(x)。
L = 1.0 # in meters
step = 1001.0 # how many steps to take (dx)
q = 0.02 # constant [N/m]
integralOfW = 0.0; # instead of w(0) enter the boundary condition value for w(0)
result = []
for i in range(int(L*step)):
x= i/step
w = (q/24.0*pow(x,4) - 0.02/12.0*pow(x,3) + 0.01/12*pow(L,3)*x)/step # current w fragment
# add up fragments of w for integral calculation
integralOfW += w
# add current value of w(x) to result list for plotting
result.append(w*step);
print("The integral of w(x) over the beam is:")
print(integralOfW)
输出：
The integral of w(x) over the beam is:
0.00016666652805511192
现在比较两种方法
射击法与数值法结果对比
w(x)在横梁上的积分：
Shooting method -> 0.000135085272117
Numerical method -> 0.00016666652805511192
这是一个很好的匹配，现在让我们看看图表：
从图中可以更明显地看出我们匹配得很好并且 shooting method 的结果是正确的。
要获得更好的 shooting method 结果，请将 xSteps 和 u2_u4_maxNumbers 增加到更大的数字，您还可以缩小 u2Numbers 和 u4Numbers到 same set size but a smaller interval（大约是上一个程序 运行 的最佳结果）。请记住，将 xSteps 和 u2_u4_maxNumbers 设置得太高会导致您的程序长时间 运行。
您需要将 ODE 转换为一阶系统，设置 u0=w 一个可能且通常使用的系统是
u0'=u1,
u1'=u2,
u2'=u3,
u3'=q(x)
这可以实现为
def ODEfunc(u,x): return [ u[1], u[2], u[3], q(x) ]
然后制作一个以实验初始条件和returns第二个边界条件的分量为准的函数
def shoot(u01, u03): return odeint(ODEfunc, [0, u01, 0, u03], [0, l])[-1,[0,2]]
现在你有一个包含两个分量的两个变量的函数，你需要用通常的方法求解这个 2x2 系统。由于系统是线性的，所以射击函数也是线性的，你只需要找到系数并求解得到的线性系统。