python code to solve the following initial value problem ordinary differential equation using Euler method over the interval (0 10)

Question

i have this question

Write a python code to solve the following initial value problem ordinary differential equation using Euler method over the interval (0 10) with 10 time steps. A) y'= -y -y^2 ; y(0)=1 If that exact solution was y(t) = 1/(-1+2e^t) What is the absolute error at y(10). now i have write this code

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())

the expected output is this

and this is the actual result that I get

i need to set the first index of y list to 1 then fill the rest of list by the for loop, how can i code this?

What is your expected output with an example input? What does your code return? — TaxpayersMoney, Jan 15 '22 at 15:47
The equations are hard to read. Do you know about LaTex? It would help. It's best to solve these in closed form if you can before you attempt the numerical solution. You should post the results you get so we can see the difference between expected and actual. — duffymo, Jan 15 '22 at 15:47
@TaxpayersMoney hi, I have attached a photo of the expected output — sajeda khalid, Jan 15 '22 at 16:43
@duffymo hi, I have attached the expected and the actual output — sajeda khalid, Jan 15 '22 at 16:47
Fix your equations. One per line. Put them in a code section. — duffymo, Jan 15 '22 at 16:48
@duffymo is it fine now ? , sorry this is my first time asking a Question here — sajeda khalid, Jan 15 '22 at 16:57
No. This makes no sense to me whatsoever: y'= -y -y?: What is up with the question and colon? Why are there two terms? Did you mean y' = -y - y^2? — duffymo, Jan 15 '22 at 16:58

John Coleman · Answer 1 · 2022-01-16T10:18:35.707

Your output is essentially correct, but shifted by 1 (so e.g. your y10 is what the intended output calls y9) and not rounded to 4 decimal places.

There are 10 updates to y but you want to print 11 values. The way to do that is to either print the first y before the loop or print the final y after the loop. The following code shows the second approach:

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}')

Note that this code simply using a scalar variable y as the variable in Euler's method (rather than an entry in an array. Ultimately a matter of taste, but the code seems cleaner that way.

The output is:

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

which matches the expected output.

I think your code is terrific, but you've only advanced the time by 10*h = 0.001 after ten steps. y(10) requires 100,000 steps of size h = 0.0001, no? You're using explicit Euler integration, which might explain why you're using such a small time step (stability). The expected result for y(10) = 0.000022700. — duffymo, Jan 18 '22 at 20:28

duffymo · Answer 2 · 2022-01-18T21:26:20.810

It helps if you can know the closed form solution before you begin.

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

This is a Bernoulli equation.

The closed form solution is:

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

Applying the initial condition to solve for the constant C1:

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

Solving for C1:

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

Substituting gives the final solution:

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

Plot the closed form solution and use that to assess your numerical integration solution.

You can get the exact solution at y(10) to compare to your numerical result:

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

I love Kotlin. I think matches Python's simplicity and succinctness but gives me more b/c of Java and the JVM. Here's what a Kotlin implementation might look like:

    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)

You can get better and better agreement with the closed form solution by taking smaller steps. Play around with the dt value to see what I mean.

You're using an explicit Euler integration scheme, which can suffer from stability limits on step sizes. You can choose more accurate and stable integration schemes like 5th order Runge Kutta.

python code to solve the following initial value problem ordinary differential equation using Euler method over the interval (0 10)

2 Answers2