Subversion Repositories wimsdev

Rev

Rev 13592 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 13592 Rev 16252
Line 26... Line 26...
26 
[1-x,x^2-y^2],0,[3,0],0.1,100,data
 26 
[1-x,x^2-y^2],0,[3,0],0.1,100,data
27 
 
 27 
 
28 
!exit
 28 
!exit
29 
:proc
 29 
:proc
30 
!distribute items $wims_read_parm into slib_rk_f,slib_rk_t0,slib_rk_y0,slib_rk_h,slib_rk_N,slib_rk_option
 30 
!distribute items $wims_read_parm into slib_rk_f,slib_rk_t0,slib_rk_y0,slib_rk_h,slib_rk_N,slib_rk_option
31 
 
 31 
!reset slib_out
32 
!if $slib_rk_h=NaN
 32 
!if $slib_rk_h=NaN
33 
  !reset slib_rk_h
 33 
  !reset slib_rk_h
34 
!endif
 34 
!endif
35 
!default slib_rk_h=0.1
 35 
!default slib_rk_h=0.1
36 
!default slib_rk_N=10
 36 
!default slib_rk_N=10
Line 54... Line 54...
54 
!for slib_rk_i=1 to $slib_rk_cnt
 54 
!for slib_rk_i=1 to $slib_rk_cnt
55 
  !if $test=1
 55 
  !if $test=1
56 
    slib_rk_tmp=$(slib_rk_y0[$slib_rk_i;])
 56 
    slib_rk_tmp=$(slib_rk_y0[$slib_rk_i;])
57 
  !else
 57 
  !else
58 
    slib_rk_tmp=[$(slib_rk_y0[$slib_rk_i;])]
 58 
    slib_rk_tmp=[$(slib_rk_y0[$slib_rk_i;])]
59 
  !endif
 59 
  !endif
60 
  !!slib_rk_program=!record 0 of gp/rungekutta.gp
 60 
  !!slib_rk_program=!record 0 of gp/rungekutta.gp
61 
 !set slib_rk_data=!exec pari rungekutta(f,t0,y0,h,N)=\
 61 
  !set slib_rk_data=!exec pari rungekutta(f,t0,y0,h,N)=\
62 
  {\
 62 
  {\
63 
  my(y,t,y1,t1);\
 63 
    my(y,t,y1,t1);\
64 
  my(X,Y);\
 64 
    my(X,Y);\
65 
  y1=y0; t1=t0;\
 65 
    y1=y0; t1=t0;\
66 
  X=vector(N); Y=vector(N);\
 66 
    X=vector(N); Y=vector(N);\
67 
  X[1]=t0; Y[1]=y0;\
 67 
    X[1]=t0; Y[1]=y0;\
68 
  for(j=2,N,\
 68 
    for(j=2,N,\
69 
      my(k1,k2,k3,k4);\
 69 
      my(k1,k2,k3,k4);\
70 
      y=y1; t=t1;\
 70 
      y=y1; t=t1;\
71 
      k1 = h*f(t, y);\
 71 
      k1 = h*f(t, y);\
72 
      k2 = h*f(t + h/2 , y + k1/2);\
 72 
      k2 = h*f(t + h/2 , y + k1/2);\
73 
      k3 = h*f(t + h/2 , y + k2/2);\
 73 
      k3 = h*f(t + h/2 , y + k2/2);\
74 
      k4 = h*f(t + h , y + k3);\
 74 
      k4 = h*f(t + h , y + k3);\
75 
      t1 = t + h;\
 75 
      t1 = t + h;\
76 
      y1 = y + 1/6 * (k1 + 2*k2 + 2*k3 + k4);\
 76 
      y1 = y + 1/6 * (k1 + 2*k2 + 2*k3 + k4);\
77 
      ;X[j]=t1;Y[j]=y1);\
 77 
      ;X[j]=t1;Y[j]=y1);\
78 
  [X,Y]\
 78 
    [X,Y]\
79 
  };\
 79 
    };\
80 
  f(t,y)={$slib_rk_f}; \
 80 
    f(t,y)={$slib_rk_f}; \
81 
  rungekutta(f,$slib_rk_t0,$slib_rk_tmp,$slib_rk_h,$slib_rk_N)
 81 
    rungekutta(f,$slib_rk_t0,$slib_rk_tmp,$slib_rk_h,$slib_rk_N)
82 
  !distribute item $slib_rk_data into slib_rk_datat, slib_rk_datay
 82 
  !distribute item $slib_rk_data into slib_rk_datat, slib_rk_datay
83 
  !set slib_rk_datat=!declosing $slib_rk_datat
 83 
  !set slib_rk_datat=!declosing $slib_rk_datat
84 
  !set slib_rk_datay=!declosing $slib_rk_datay
 84 
  !set slib_rk_datay=!declosing $slib_rk_datay
85 
  !reset slib_rk_d
 85 
  !reset slib_rk_d
86 
  !if $test=1
 86 
  !if $test=1
Line 107... Line 107...
107 
!endif
 107 
!endif
108 
 
 108 
 
109 
!set slib_rk_draw=xrange $(slib_rk_rangex[1])*1.1,$(slib_rk_rangex[2])*1.1\
 109 
!set slib_rk_draw=xrange $(slib_rk_rangex[1])*1.1,$(slib_rk_rangex[2])*1.1\
110 
 yrange $(slib_rk_rangey[1])*1.1, $(slib_rk_rangey[2])*1.1
 110 
 yrange $(slib_rk_rangey[1])*1.1, $(slib_rk_rangey[2])*1.1
111 
!for slib_rk_i=1 to $slib_rk_cnt
 111 
!for slib_rk_i=1 to $slib_rk_cnt
- 
 
 112 
  !set slib_tmp=!nonempty lines $(slib_out[$slib_rk_i;])
112 
    slib_rk_draw=!append line polyline black,$(slib_out[$slib_rk_i;])\
 113 
    slib_rk_draw=!append line polyline black,$slib_tmp\
113 
    fcircle $(slib_out[$slib_rk_i;1]),$(slib_out[$slib_rk_i;2]),10,blue to $slib_rk_draw
 114 
    fcircle $(slib_out[$slib_rk_i;1]),$(slib_out[$slib_rk_i;2]),10,blue to $slib_rk_draw
114 
!next
 115 
!next
115 
slib_rk_draw=$slib_rk_draw\
 116 
slib_rk_draw=$slib_rk_draw\
116 
 arrow 0,$(slib_rk_rangey[1]),0,$(slib_rk_rangey[2]),6,black\
 117 
 arrow 0,$(slib_rk_rangey[1]),0,$(slib_rk_rangey[2]),6,black\
117 
 arrow $(slib_rk_rangex[1]),0,$(slib_rk_rangex[2]),0,6,black\
 118 
 arrow $(slib_rk_rangex[1]),0,$(slib_rk_rangex[2]),0,6,black\
118 
 
 - 
 
119 
!if url iswordof $slib_rk_option or html iswordof $slib_rk_option
 119 
!if url iswordof $slib_rk_option or html iswordof $slib_rk_option
120 
  !insdraw $slib_rk_draw
 120 
  !insdraw $slib_rk_draw
121 
  slib_out=$ins_url
 121 
  slib_out=$ins_url
122 
  !if html iswordof $slib_rk_option
 122 
  !if html iswordof $slib_rk_option
123 
    slib_out=<img src="$ins_url" alt="" />
 123 
    slib_out=<img src="$ins_url" alt="" />