Question: How do i resolve my code which encountering error in a loop for repeated differentiation?

i am trying to make code for automatic differentiation method. kindly help me out with that. i have asked chatgpt too. i am going to attach my file alongwith chatgpt code kindly help meautomatic_differentiation.mw
 

restart

v := 1; a := 2; t := 0.1e-2; dt := 0.1e-3; N := 40; h := 1/40; K := 4

NULL

NULL

for i from 0 to N do x[i] := i*h end do

NULL

initial_conditions := []; for i to N do initial_conditions := [op(initial_conditions), 2*v*Pi*sin(Pi*x[i])/(a+cos(Pi*x[i]))] end do

u := proc (i) local u_x, u_xx, expr, j; u_x := (1/2)*(u[i+1]-u[i-1])/h; u_xx := (u[i-1]-2*u[i]+u[i+1])/h^2; expr := -alpha*u[i]*u_x+v*u_xx; expr end proc
NULL

NULL

odes := [seq(u(i, [seq(u[j], j = 1 .. N-1)]), i = 1 .. N-1)]

for i to N-1 do assign(o[i] = odes[i]) end do

for i to N-1 do printf("u_%d = %s\n", i, convert(u(i), string)) end do

u_1 = -alpha*u[1]*(20*u[2]-20*u[0])+1600*u[0]-3200*u[1]+1600*u[2]
u_2 = -alpha*u[2]*(20*u[3]-20*u[1])+1600*u[1]-3200*u[2]+1600*u[3]
u_3 = -alpha*u[3]*(20*u[4]-20*u[2])+1600*u[2]-3200*u[3]+1600*u[4]
u_4 = -alpha*u[4]*(20*u[5]-20*u[3])+1600*u[3]-3200*u[4]+1600*u[5]
u_5 = -alpha*u[5]*(20*u[6]-20*u[4])+1600*u[4]-3200*u[5]+1600*u[6]
u_6 = -alpha*u[6]*(20*u[7]-20*u[5])+1600*u[5]-3200*u[6]+1600*u[7]
u_7 = -alpha*u[7]*(20*u[8]-20*u[6])+1600*u[6]-3200*u[7]+1600*u[8]
u_8 = -alpha*u[8]*(20*u[9]-20*u[7])+1600*u[7]-3200*u[8]+1600*u[9]
u_9 = -alpha*u[9]*(20*u[10]-20*u[8])+1600*u[8]-3200*u[9]+1600*u[10]
u_10 = -alpha*u[10]*(20*u[11]-20*u[9])+1600*u[9]-3200*u[10]+1600*u[11]
u_11 = -alpha*u[11]*(20*u[12]-20*u[10])+1600*u[10]-3200*u[11]+1600*u[12]
u_12 = -alpha*u[12]*(20*u[13]-20*u[11])+1600*u[11]-3200*u[12]+1600*u[13]
u_13 = -alpha*u[13]*(20*u[14]-20*u[12])+1600*u[12]-3200*u[13]+1600*u[14]
u_14 = -alpha*u[14]*(20*u[15]-20*u[13])+1600*u[13]-3200*u[14]+1600*u[15]
u_15 = -alpha*u[15]*(20*u[16]-20*u[14])+1600*u[14]-3200*u[15]+1600*u[16]
u_16 = -alpha*u[16]*(20*u[17]-20*u[15])+1600*u[15]-3200*u[16]+1600*u[17]
u_17 = -alpha*u[17]*(20*u[18]-20*u[16])+1600*u[16]-3200*u[17]+1600*u[18]
u_18 = -alpha*u[18]*(20*u[19]-20*u[17])+1600*u[17]-3200*u[18]+1600*u[19]
u_19 = -alpha*u[19]*(20*u[20]-20*u[18])+1600*u[18]-3200*u[19]+1600*u[20]
u_20 = -alpha*u[20]*(20*u[21]-20*u[19])+1600*u[19]-3200*u[20]+1600*u[21]
u_21 = -alpha*u[21]*(20*u[22]-20*u[20])+1600*u[20]-3200*u[21]+1600*u[22]
u_22 = -alpha*u[22]*(20*u[23]-20*u[21])+1600*u[21]-3200*u[22]+1600*u[23]
u_23 = -alpha*u[23]*(20*u[24]-20*u[22])+1600*u[22]-3200*u[23]+1600*u[24]
u_24 = -alpha*u[24]*(20*u[25]-20*u[23])+1600*u[23]-3200*u[24]+1600*u[25]
u_25 = -alpha*u[25]*(20*u[26]-20*u[24])+1600*u[24]-3200*u[25]+1600*u[26]
u_26 = -alpha*u[26]*(20*u[27]-20*u[25])+1600*u[25]-3200*u[26]+1600*u[27]
u_27 = -alpha*u[27]*(20*u[28]-20*u[26])+1600*u[26]-3200*u[27]+1600*u[28]
u_28 = -alpha*u[28]*(20*u[29]-20*u[27])+1600*u[27]-3200*u[28]+1600*u[29]
u_29 = -alpha*u[29]*(20*u[30]-20*u[28])+1600*u[28]-3200*u[29]+1600*u[30]
u_30 = -alpha*u[30]*(20*u[31]-20*u[29])+1600*u[29]-3200*u[30]+1600*u[31]
u_31 = -alpha*u[31]*(20*u[32]-20*u[30])+1600*u[30]-3200*u[31]+1600*u[32]
u_32 = -alpha*u[32]*(20*u[33]-20*u[31])+1600*u[31]-3200*u[32]+1600*u[33]
u_33 = -alpha*u[33]*(20*u[34]-20*u[32])+1600*u[32]-3200*u[33]+1600*u[34]
u_34 = -alpha*u[34]*(20*u[35]-20*u[33])+1600*u[33]-3200*u[34]+1600*u[35]
u_35 = -alpha*u[35]*(20*u[36]-20*u[34])+1600*u[34]-3200*u[35]+1600*u[36]
u_36 = -alpha*u[36]*(20*u[37]-20*u[35])+1600*u[35]-3200*u[36]+1600*u[37]
u_37 = -alpha*u[37]*(20*u[38]-20*u[36])+1600*u[36]-3200*u[37]+1600*u[38]
u_38 = -alpha*u[38]*(20*u[39]-20*u[37])+1600*u[37]-3200*u[38]+1600*u[39]
u_39 = -alpha*u[39]*(20*u[40]-20*u[38])+1600*u[38]-3200*u[39]+1600*u[40]

 

eval_derivatives := proc (k, u) local T1i, T2i, T3i, T1ik, T2ik, T3ik, ui_k, i, j; T1i := Array(0 .. N); T2i := Array(0 .. N); T3i := Array(0 .. N); ui_k := Array(0 .. N); for i to N-1 do T1i[i] := (1/2)*(u[i+1]-u[i-1])/h; T2i[i] := u[i]*T1i[i]; T3i[i] := (u[i-1]-2*u[i]+u[i+1])/h^2; ui_k[i] := v*T3i[i]-T2i[i] end do; T1ik := Array(0 .. N); T2ik := Array(0 .. N); T3ik := Array(0 .. N); for i to N-1 do T1ik[i] := (1/2)*(u[i+1, k]-u[i-1, k])/h; T2ik[i] := add(u[i, j]*T1i[i, k-j], j = 0 .. k); T3ik[i] := (u[i-1, k]-2*u[i, k]+u[i+1, k])/h^2; ui_k[i] := (v*T3ik[i]-T2ik[i])/(k+1) end do; return ui_k end proc
``

proc (k, u) local T1i, T2i, T3i, T1ik, T2ik, T3ik, ui_k, i, j; T1i := Array(0 .. N); T2i := Array(0 .. N); T3i := Array(0 .. N); ui_k := Array(0 .. N); for i to N-1 do T1i[i] := (1/2)*(u[i+1]-u[i-1])/h; T2i[i] := u[i]*T1i[i]; T3i[i] := (u[i-1]-2*u[i]+u[i+1])/h^2; ui_k[i] := v*T3i[i]-T2i[i] end do; T1ik := Array(0 .. N); T2ik := Array(0 .. N); T3ik := Array(0 .. N); for i to N-1 do T1ik[i] := (1/2)*(u[i+1, k]-u[i-1, k])/h; T2ik[i] := add(u[i, j]*T1i[i, k-j], j = 0 .. k); T3ik[i] := (u[i-1, k]-2*u[i, k]+u[i+1, k])/h^2; ui_k[i] := (v*T3ik[i]-T2ik[i])/(k+1) end do; return ui_k end proc

(1)

u := Array(1 .. N, 0 .. K); for i to N do u[i, 0] := evalf(initial_conditions[i]) end do

Array(%id = 36893489585014251020)

 

.1644933719

 

.3289856976

 

.4934717144

 

.6579375992

 

.8223563570

 

.9866828090

 

1.150848028

 

1.314753051

 

1.478261707

 

1.641192349

 

1.803308276

 

1.964306617

 

2.123805434

 

2.281328760

 

2.436289364

 

2.587968970

 

2.735495794

 

2.877819368

 

3.013682762

 

3.141592654

 

3.259788188

 

3.366210070

 

3.458472370

 

3.533840560

 

3.589220824

 

3.621167336

 

3.625916040

 

3.599455182

 

3.537643690

 

3.436388151

 

3.291886154

 

3.100937330

 

2.861312732

 

2.572157998

 

2.234388242

 

1.851015873

 

1.427342882

 

.9709526944

 

.4914580366

 

0.

(2)

for k to K do ui_k := eval_derivatives(k, u); for i to N-1 do u[i, k] := ui_k[i] end do end do

Error, (in eval_derivatives) Array index out of range

 

NULL

"     local u_new, i, k;      u_new := Array(1..N);      for i from 1 to N do          u_new[i] := evalf(add(u[i, k] * dt^k / factorial(k), k = 0 .. K));      end do;      return u_new;  end proc;"

Error, unable to parse

"     local u_new, i, k;   u_new := Array(1..N);      for i from 1 to N do   u_new[i] := evalf(add(u[i, k] * dt^k / factorial(k), k = 0 .. K));      end do;   return u_new;  end proc;"

 

num_steps := round(0.1e-2/dt); for step to num_steps do u_new := advance_solution(t, dt, K, u); for i to N do u[i, 0] := u_new[i] end do; for k to K do ui_k := eval_derivatives(k, u); for i to N do u[i, k] := ui_k[i] end do end do; t := t+dt end do; print(u_new)

num_steps := 10

 

advance_solution(0.1e-2, 0.1e-3, 4, Array(%id = 36893489585014251020))

(3)

NULL


 

Download automatic_differentiation.mw

resolve the issue

Please Wait...