stabilisation par RK4-Adams4

176 days ago by mhbahichem

import numpy as np def f(t,y): res=[y[1],-(y[0]^2)] return(np.array(res)) 
       
def RK4(a,b,N,ya): y=np.empty([2,N+1]) #tableau 2 lignes N+1 colonnes h=(b-a)/N y[:,0]=ya for i in range(N): t=a+i*h tdemi=t+h/2 r1=f(t,y[:,i]) r2=f(tdemi,y[:,i]+(h/2)*r1) r3=f(tdemi,y[:,i]+(h/2)*r2) r4=f(t,y[:,i]+h*r3) y[:,i+1]=y[:,i]+(h/6)*(r1+2*r2+2*r3+r4) return (y) 
       
RK4(0,5,10,1) 
       
array([[  1.00000000e+00,   1.33447266e+00,   1.24812557e+00,
          7.90350203e-01,   1.68766077e-01,  -4.76809912e-01,
         -1.19828700e+00,  -2.32649094e+00,  -5.04840935e+00,
         -1.65036736e+01,  -1.91557328e+02],
       [  1.00000000e+00,   2.79373169e-01,  -6.01953019e-01,
         -1.15220442e+00,  -1.28485955e+00,  -1.31399740e+00,
         -1.67595853e+00,  -3.18259805e+00,  -9.45245958e+00,
         -5.73215667e+01,  -2.86171751e+03]])
array([[  1.00000000e+00,   1.33447266e+00,   1.24812557e+00,
          7.90350203e-01,   1.68766077e-01,  -4.76809912e-01,
         -1.19828700e+00,  -2.32649094e+00,  -5.04840935e+00,
         -1.65036736e+01,  -1.91557328e+02],
       [  1.00000000e+00,   2.79373169e-01,  -6.01953019e-01,
         -1.15220442e+00,  -1.28485955e+00,  -1.31399740e+00,
         -1.67595853e+00,  -3.18259805e+00,  -9.45245958e+00,
         -5.73215667e+01,  -2.86171751e+03]])
def Adams(a,b,N,ya): y=np.empty([2,N+1]) h=(b-a)/N y[:,0]=ya rk4=RK4(a,a+2*h,2,ya) y[:,0:3]=rk4 for i in range(2,N): t=a+i*h y[:,i+1]=y[:,i]+(h/24)*(9*f(t,y[:,i+1])+19*f(t,y[:,i])-5*f(t,y[:,i-1])+f(t,y[:,i-2])) return(y) 
       
Adams(0,5,50,[1]) 
       
Traceback (click to the left of this block for traceback)
...
TypeError: only length-1 arrays can be converted to Python scalars
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_137.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("bT1BZGFtcygwLDUsNTAsWzFdKQpsaXN0X3Bsb3QobGlzdChtKSk="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmp_k4xKz/___code___.py", line 4, in <module>
    exec compile(u'list_plot(list(m))
  File "", line 1, in <module>
    
  File "/sagenb/sage_install/sage-4.7.2/local/lib/python2.6/site-packages/sage/misc/decorators.py", line 534, in wrapper
    return func(*args, **options)
  File "/sagenb/sage_install/sage-4.7.2/local/lib/python2.6/site-packages/sage/plot/plot.py", line 3656, in list_plot
    P = point(data, **kwargs)
  File "/sagenb/sage_install/sage-4.7.2/local/lib/python2.6/site-packages/sage/plot/point.py", line 312, in point
    return point3d(points, **kwds)
  File "/sagenb/sage_install/sage-4.7.2/local/lib/python2.6/site-packages/sage/plot/plot3d/shapes2.py", line 1022, in point3d
    A = sum([Point(z, size, **kwds) for z in v])
  File "/sagenb/sage_install/sage-4.7.2/local/lib/python2.6/site-packages/sage/plot/plot3d/shapes2.py", line 680, in __init__
    self.loc = (float(center[0]), float(center[1]), float(center[2]))
TypeError: only length-1 arrays can be converted to Python scalars