الحل العددي للنماذج الرياضية للكائنات المعطاة بواسطة أنظمة المعادلات التفاضلية

from numpy import* import matplotlib.pyplot as plt from scipy.integrate import * def odein(): #dy1/dt=y2 #dy2/dt=y1**2+1: def f(y,t): return y**2+1 t =arange(0,1,0.01) y0 =0.0 y=odeint(f, y0,t) y = array(y).flatten() return y,t def oden(): f = lambda t, y: y**2+1 ODE=ode(f) ODE.set_integrator('dopri5') ODE.set_initial_value(0, 0) t=arange(0,1,0.01) z=[] t=arange(0,1,0.01) for i in arange(0,1,0.01): ODE.integrate(i) q=ODE.y z.append(q[0]) return z,t def rungeKutta(f, to, yo, tEnd, tau): def increment(f, t, y, tau): if z==1: k0 =tau* f(t,y) k1 =tau* f(t+tau/2.,y+k0/2.) k2 =tau* f(t+tau/2.,y+k1/2.) k3 =tau* f(t+tau, y + k2) return (k0 + 2.*k1 + 2.*k2 + k3) / 6. elif z==0: k1=tau*f(t,y) k2=tau*f(t+(1/4)*tau,y+(1/4)*k1) k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2) k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3) k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4) k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5) return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6 t = [] y= [] t.append(to) y.append(yo) while to < tEnd: tau = min(tau, tEnd - to) yo = yo + increment(f, to, yo, tau) to = to + tau t.append(to) y.append(yo) return array(t), array(y) def f(t, y): f = zeros([1]) f[0] = y[0]**2+1 return f to = 0. tEnd = 1 yo = array([0.]) tau = 0.01 z=1 t, yn = rungeKutta(f, to, yo, tEnd, tau) y1n=[i[0] for i in yn] plt.figure() plt.title("   (..- u(t)=tan(t)) \n\ du/dt=u**2+1 cu(0)=0  t>0") plt.plot(t,abs(array(y1n)-array(tan(t))),label=' — \n\   -   ') plt.xlabel('') plt.ylabel(' .') plt.legend(loc='best') plt.grid(True) z=0 t, ym = rungeKutta(f, to, yo, tEnd, tau) y1m=[i[0] for i in ym] plt.figure() plt.title("   (..- u(t)=tan(t)) \n\ du/dt=u**2+1 cu(0)=0  t>0") plt.plot(t,abs(array(y1m)-array(tan(t))),label=' ——  \n\   -   ') plt.xlabel('') plt.ylabel(' .') plt.legend(loc='best') plt.grid(True) plt.figure() plt.title("    (..- u(t)=tan(t)) \n\ du/dt=u**2+1 cu(0)=0  t>0") y,t=odein() plt.plot(t,abs(array(tan(t))-array(y)),label=' odein') plt.xlabel('') plt.ylabel(' .') plt.legend(loc='best') plt.grid(True) plt.figure() plt.title("    (..- u(t)=tan(t)) \n\ du/dt=u**2+1 cu(0)=0  t>0") z,t=oden() plt.plot(t,abs(tan(t)-z),label=' ode  ——  \n\  ') plt.xlabel('') plt.ylabel(' .') plt.legend(loc='best') plt.grid(True) plt.show() 

 import numpy as np import matplotlib.pyplot as plt import time start = time.time() from scipy.integrate import ode ts = [ ] ys = [ ] FlightTime, Distance, Height =0,0,0 y4old=0 def fout(t, y):#   global FlightTime, Distance, Height,y4old ts.append(t) ys.append(list(y.copy())) y1, y2, y3, y4 = y if y4*y4old<=0: #    Height=y3 if y4<0 and y3<=0.0: #    FlightTime=t Distance=y1 return -1 y4old=y4 #      def f(t, y, k): #    k g=9.81 y1, y2, y3, y4 = y return [y2,-k*y2*np.sqrt(y2**2+y4**2), y4,-k*y4*np.sqrt(y2**2+y4**2)-g] tmax=1.41 #     alph=np.pi/4 #    v0=10.0 #   K=[0.1,0.2,0.3,0.5] #    y0,t0=[0, v0*np.cos(alph), 0, v0*np.sin(alph)], 0 #   ODE=ode(f) ODE.set_integrator('dopri5', max_step=0.01) ODE.set_solout(fout) fig, ax = plt.subplots() fig.set_facecolor('white') for k in K: #     ts, ys = [ ],[ ] ODE.set_initial_value(y0, t0) #    ODE.set_f_params(k) #    k #   f(t,y,k)     ODE.integrate(tmax) #   print('Flight time = %.4f Distance = %.4f Height =%.4f '% (FlightTime,Distance,Height)) Y=np.array(ys) plt.plot(Y[:,0],Y[:,2],linewidth=3,label='k=%.1f'% k) stop = time.time() plt.title("      \n    ode   dopri5 ") print ("   : %f"%(stop-start)) plt.grid(True) plt.xlim(0,8) plt.ylim(-0.1,2) plt.legend(loc='best') plt.show() 

 from numpy import* import matplotlib.pyplot as plt import time start = time.time() def rungeKutta(f, to, yo, tEnd, tau): def increment(f, t, y, tau):#     ——. k1=tau*f(t,y) k2=tau*f(t+(1/4)*tau,y+(1/4)*k1) k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2) k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3) k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4) k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5) return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6 t = []#   t y= []#   y t.append(to)#   t   to y.append(yo)#   y   yo while to < tEnd:#     t,y tau = min(tau, tEnd - to)#   tau yo = yo + increment(f, to, yo, tau) #     t0,y0    to = to + tau #   t.append(to) #   t y.append(yo) #   y return array(t), array(y) def f(t, y): #      f = zeros([4]) f[0]=y[1] f[1]=-k*y[1]*sqrt(y[1]**2+y[3]**2) f[2]=y[3] f[3]=-k*y[3]*sqrt(y[1]**2+y[3]**2) -g if y[3]<0 and y[2]<=0.0: #    return -1 return f to = 0#     tEnd = 1.41#     alph=pi/4#    v0=10.0 #   K=[0.1,0.2,0.3,0.5]#     g=9.81 yo = array([0.,v0*cos(alph),0.,v0*sin(alph)]) #   tau =0.01#  for i in K: #      k=i t, y = rungeKutta(f, to, yo, tEnd, tau) y1=array([i[0] for i in y]) #     y y3=array([i[2] for i in y]) #    ""     s,h,t plt.plot(y1,y3,linewidth=2,label='k=%.1f h=%.3f s=%.2f t=%s' % (k,max(y3),max(y1),round(t[list(y1).index(max(y1))],3))) stop = time.time() plt.title("      \n     Python\n  —— ") print ("   : %f"%(stop-start)) plt.xlabel(' h') plt.ylabel(' s') plt.legend(loc='best') plt.xlim(0,8) plt.ylim(-0.1,2) plt.grid(True) plt.show() 

  #   from numpy import* import matplotlib.pyplot as plt import matplotlib.font_manager as fm,os import matplotlib.patches as mpatches import matplotlib.lines as mlines from scipy.integrate import odeint from scipy import linalg import time start = time.time() c1 = 1.0 c2 = 0.8 c3 = 0.5 a =0.0 b = 1.0 nn =100 initial_state_0 =array( [a, c1, 0.0, 0.0]) initial_state_I =array( [a, 0.0, 1.0, 0.0]) initial_state_II =array( [a, 0.0, 0.0, 1.0]) to = a tEnd =b N = int(nn) tau=(ba)/N def rungeKutta(f, to, yo, tEnd, tau): def increment(f, t, y, tau): k1=tau*f(t,y) k2=tau*f(t+(1/4)*tau,y+(1/4)*k1) k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2) k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3) k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4) k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5) return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6 t = [] y= [] t.append(to) y.append(yo) while to < tEnd: tau = min(tau, tEnd - to) yo = yo + increment(f, to, yo, tau) to = to + tau t.append(to) y.append(yo) return array(t), array(y) def f(t, y): global theta f = zeros([4]) f[0] = 1 f[1] = -y [1]-y[2] +theta* sin(y[0]) f[2] = -y[2]+y[3] f[3] = -y[2] return f #    -- theta = 1 theta = 1.0 yo =initial_state_0 t, y = rungeKutta(f, to, yo, tEnd, tau) y2=[i[2] for i in y] y3=[i[3] for i in y] #       # Y20 = Y2(b), Y30 = Y3(b) Y20 = y2[N-1] Y30 = y3[N-1] #    -- theta = 0,  I theta = 0.0 yo= initial_state_I t, y = rungeKutta(f, to, yo, tEnd, tau) y2=[i[2] for i in y] y3=[i[3] for i in y] #       # Y21= Y2(b), Y31 = Y3(b) Y21= y2[N-1] Y31 = y3[N-1] #    -- theta = 0,  II theta = 0.0 yo =initial_state_II t, y = rungeKutta(f, to, yo, tEnd, tau) y2=[i[2] for i in y] y3=[i[3] for i in y] #       # Y211= Y2(b), Y311 = Y3(b) Y211= y2[N-1] Y311 = y3[N-1] #    #     p2, p3 b1 = c2 - Y20 b2 = c3 - Y30 A = array([[Y21, Y211], [Y31, Y311]]) bb = array([[b1], [b2]]) #   p2, p3 = linalg.solve(A, bb) #    #  , theta = 1 theta = 1.0 yo = array([a, c1, p2, p3]) t, y = rungeKutta(f, to, yo, tEnd, tau) y0=[i[0] for i in y] y1=[i[1] for i in y] y2=[i[2] for i in y] y3=[i[3] for i in y] #  print('y0[0]=', y0[0]) print('y1[0]=', y1[0]) print('y2[0]=', y2[0]) print('y3[0]=', y3[0]) print('y0[N-1]=', y0[N-1]) print('y1[N-1]=', y1[N-1]) print('y2[N-1]=', y2[N-1]) print('y3[N-1]=', y3[N-1]) j = N xx = y0[:j] yy1 = y1[:j] yy2 = y2[:j] yy3 = y3[:j] stop = time.time() print ("   : %f"%(stop-start)) plt.subplot(2, 1, 1) plt.plot([a], [c1], 'ro') plt.plot([b], [c2], 'go') plt.plot([b], [c3], 'bo') plt.plot(xx, yy1, color='r') #  plt.plot(xx, yy2, color='g') #  plt.plot(xx, yy3, color='b') #  plt.xlabel(r'$x$') #   x   TeX plt.ylabel(r'$y_k(x)$') #   y   TeX plt.title(r'  ', color='blue') plt.grid(True) # patch_y1 = mpatches.Patch(color='red', label='$y_1$') patch_y2 = mpatches.Patch(color='green', label='$y_2$') patch_y3 = mpatches.Patch(color='blue', label='$y_3$') plt.legend(handles=[patch_y1, patch_y2, patch_y3]) ymin, ymax = plt.ylim() xmin, xmax = plt.xlim() plt.subplot(2, 1, 2) font = {'family': 'serif', 'color': 'blue', 'weight': 'normal', 'size': 12, } plt.text(0.2, 0.8, r'$\frac{dy_1}{dx}= - y_1 - y_2 + \sin(x),$', fontdict=font) plt.text(0.2, 0.6,r'$\frac{dy_2}{dx}= - y_1 + y_3,$', fontdict=font) plt.text(0.2, 0.4, r'$\frac{dy_3}{dx}= - y_2 - y_2,$', fontdict=font) plt.text(0.2, 0.2, r'$y_1(a)=c_1, ' r'\quad y_2(b)=c_2, \quad y_3(b)=c_3.$', fontdict=font) plt.subplots_adjust(left=0.15) plt.show()