问python代码在区间(0 ~ 10)上用Euler方法求解以下常微分方程初值问题
EN

Stack Overflow用户

提问于 2022-01-15 07:42:31

回答 2查看 665关注 0票数 0

我有个问题

编写一个python代码，用Euler方法在区间(0 ~ 10)上用10个时间步长求解以下初值问题--常微分方程。( A) y'= -y -y^2；y(0)=1，如果这个精确解是y(t) = 1/(-1+2e^t)，那么y(10)的绝对误差是多少。现在我已经写了这段代码

def pow_t(x,b):  
    t=1; 
    while (b):
        t*=x 
        b=b-1
    return t


def  absolute(): 
    y=[x for x in range(1,11)]
    
    h=0.0001 
    for i in range(1,10) :
       y[i]=y[i-1]+(h*(-1*y[i-1]-pow_t(y[i-1],2)))
       print("y",i,"=",y[i]) 
    
    exact = 0.0000227
    approx = y[9]
    absolute = exact - approx 
    print("abbsolute erroe = exact - approx ")
    print("abbsolute erroe = ",absolute)
   

print(absolute())

预期的输出如下

这是我得到的实际结果

我需要将y列表的第一个索引设置为1，然后通过for循环填充列表的其余部分，我如何编码呢？

python

math

differential-equations

runge-kutta

回答 2

Stack Overflow用户

发布于 2022-01-15 11:07:39

您的输出基本上是正确的，但是移动了1(例如，您的y10是预期输出所称的y9)，而不是四舍五入到小数点后4位。

y有10个更新，但要打印11个值。这样做的方法是要么在循环之前打印第一个y，要么在循环之后打印最终的y。下面的代码显示了第二种方法：

y = 1
t = 0
h = 0.0001
iterates = [1]

for i in range(10):
    print(f'y{i} = {y:.4f}')
    y = y+h*(-y-y**2)
    iterates.append(y)
    t += h
print(f'y10 = {y:.4f}')

注意，这段代码只使用标量变量y作为欧拉方法中的变量(而不是数组中的条目)。说到底，这是一个品味的问题，但这样的代码似乎更干净。

产出如下：

y0 = 1.0000
y1 = 0.9998
y2 = 0.9996
y3 = 0.9994
y4 = 0.9992
y5 = 0.9990
y6 = 0.9988
y7 = 0.9986
y8 = 0.9984
y9 = 0.9982
y10 = 0.9980

它与预期的输出相匹配。

票数 1

Stack Overflow用户

发布于 2022-01-15 09:48:14

如果您可以在开始之前知道封闭表单解决方案，这会有所帮助。

Equation: y' + y + y^2 = 0
Initial condition: y(0) = 1

这是伯努利方程。

闭形解是：

y(t) = -e^C1 / (e^C1 - e^t)

应用常量C1解的初始条件：

y(0) = -e^C1 / (e^C1 - 1) = 1

C1的解决方案：

e^C1 = 1/2
C1 = ln(1/2)

替换给出了最终的解决方案：

y(t) = 1/(2*(e^t - 1/2))

绘制封闭形式的解决方案，并使用它来评估数值积分解决方案。

您可以在y(10)得到精确的解，以便与数值结果进行比较：

y(0)  = 1
y(1)  = 0.225399674
y(2)  = 0.072578883
y(3)  = 0.025529042
y(4)  = 0.009242460
y(5)  = 0.003380362
y(6)  = 0.001240914
y(7)  = 0.000456149
y(8)  = 0.000671038
y(9)  = 0.000061709
y(10) = 1/(2*(e^(10)) - 1/2)) = 0.000022700

我爱Kotlin。我认为Python的简单性和简洁性相匹配，但给了我更多的Java和JVM的b/c。下面是Kotlin实现的样子：

    fun main() {
        var y = 1.0
        var t = 0.0
        val dt = 0.005
        val tmax = 10.0
        do {
            println("y(${t.format(5)}) = ${y.format(10)}")
            y += dt*(-y-y*y)
            t += dt
        } while (t < tmax)
    }

fun Double.format(digits: Int) = "%.${digits}f".format(this)