Advanced Search

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

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

Khair Muhammad Saraz 40 Maple
differentiation
August 07 2024
0

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...