• Explore
• My scripts
• Calculator

Schéma d’Euler, schéma de Runge-Kutta (RK1=Euler, RK2 = point milieu, RK4) Exemple y’=sin(xy) avec (x0,y0)=(1,2), valeur en x=3 en 1000 pas s’obtient par RK2(lambda x,y: sin(xy),1,2,3,1000) https://fr.wikipedia.org/wiki/M%C3%A9thodes_de_Runge-Kutta

from math import *
from matplotlib.pyplot import *

def Euler(f, x0, y0, xn, n=100):
  x = x0
  y = y0
  h=(xn-x0)/n
  for _i in range(n):
    k1 = h * f(x, y)
    y = y + k1
    x = x + h
  return y

def RK2(f, x0, y0, xn, n=100):
  x = x0
  y = y0
  h=(xn-x0)/n
  for _i in range(n):
    k1 = h * f(x, y)
    k2 = h * f(x + h/2, y + k1/2)
    y = y + k2
    x = x + h
  return y

def RK4(f, x0, y0, xn, n=50):
  xs,ys=[],[] # pour tracer
  my,My=y0,y0 
  x = x0
  y = y0
  h=(xn-x0)/n
  for _i in range(n):
    xs.append(x)
    ys.append(y)
    if(my>y): my=y
    if(My<y): My=y
    k1 = h * f(x, y)
    k2 = h * f(x + h/2, y + k1/2)
    k3 = h * f(x + h/2, y + k2/2)
    k4 = h * f(x + h, y + k3)
    y = y + (k1+2*k2+2*k3+k4)/6
    x = x + h
  plot(xs,ys)
  champ(f,x0, xn, my, My)
  show()
  return y

def champ(f,x0,x1,y0,y1,n=10):
  axis((x0, x1, y0, y1))
  grid()
  hx=(x1-x0)/n
  hy=(y1-y0)/n
  r2=(hy/hx)**2
  for i in range(n):
    for j in range (n):
      x = x0+i*hx
      y = y0+j*hy
      fxy=f(x,y)
      l=1/sqrt(r2+fxy**2)/2
      plot((x,x+hy*l),(y,y+hy*fxy*l))

RK4(lambda x, y: -2*x*y, -2, 0.1, 2)