Advanced Search

Question: differential quadrature method stuck at rk4 method coding after successfully derived system of odes using dqm

Question:differential quadrature method stuck at rk4 method coding after successfully derived system of odes using dqm

Khair Muhammad Saraz 35 Maple 2022
procedure loop homework numeric differential-equations plotting
May 20 2024
1

im trying to write differential quadrature method i have successfully  generated system of equations using code but now im stuck with rk4 method to solve these equation in a loop kindly give me suggestions or modify the code

> 

restart
> 

with(orthopoly)
> 

N := 5

5

 (1)
> 

P := proc (n, x) options operator, arrow; orthopoly[P](n, 2*x-1) end proc
> 

P_N := simplify(P(N, x))

252*x^5-630*x^4+560*x^3-210*x^2+30*x-1

 (2)
> 

localroots := fsolve(P_N = 0, x, complex)

HFloat(0.04691007703066799), HFloat(0.23076534494715836), HFloat(0.4999999999999997), HFloat(0.7692346550528439), HFloat(0.9530899229693298)

 (3)
> 

for i to N do assign(x[i] = localroots[i]) end do
> 

PRime := diff(P_N, x)

1260*x^4-2520*x^3+1680*x^2-420*x+30

 (4)
> 

for i to N do for j to N do if i <> j then a[i, j] := (eval(PRime, x = x[i]))/((x[i]-x[j])*(eval(PRime, x = x[j]))) else a[i, j] := (1-2*x[i])/(2*x[i]*(x[i]-1)) end if end do end do
> 

u := proc (i) local u_x, expr, j; u_x := add(a[i, j]*u[j], j = 1 .. N); expr := -u[i]*u_x+x[i]; expr end proc; odes := [seq(u(i, [seq(u[j], j = 1 .. N)]), i = 1 .. N)]; for i to N do assign(o[i] = odes[i]) end do; for i to N do printf("u_%d = %s\n", i, convert(u(i), string)) end do

u_1 = -u[1]*(-10.1340811913091*u[1]+15.4039041703445*u[2]-8.08708750877537*u[3]+3.92079823166652*u[4]-1.10353370192667*u[5])+.046910077030668
u_2 = -u[2]*(-1.92051204726391*u[1]-1.51670643433575*u[2]+4.80550130432862*u[3]-1.85711605328765*u[4]+.488833230558754*u[5])+.230765344947158
u_3 = -u[3]*(.602336319455661*u[1]-2.87077648466948*u[2]-1.11022302462516e-015*u[3]+2.87077648466942*u[4]-.602336319455684*u[5])+.5
u_4 = -u[4]*(-.488833230558737*u[1]+1.85711605328767*u[2]-4.8055013043286*u[3]+1.51670643433578*u[4]+1.92051204726404*u[5])+.769234655052844
u_5 = -u[5]*(1.1035337019266*u[1]-3.92079823166643*u[2]+8.08708750877511*u[3]-15.4039041703442*u[4]+10.1340811913086*u[5])+.95308992296933

  
> 

initial_conditions := [0.1e-1, 0.2e-1, 0.3e-1, 0.4e-1, 0.5e-1]
> 

for i to N do printf("Initial condition for u_%d: u_%d(0) = %f\n", i, i, initial_conditions[i]) end do

Initial condition for u_1: u_1(0) = 0.010000
Initial condition for u_2: u_2(0) = 0.020000
Initial condition for u_3: u_3(0) = 0.030000
Initial condition for u_4: u_4(0) = 0.040000
Initial condition for u_5: u_5(0) = 0.050000

  
> 

dt := .1
> 

tf := 1
> 

"local o_old; o_old := Vector(n); #` Loop to initialize o[i]old values` for i from 1 to n do o_old[i] := initial_conditions[i]; printf("Initial condition for u_%d: u_%d(0) = %f\\n", i, i, initial_conditions[i]); end do: time := t0; while time < tf do for i from 1 to n do k[1, i] := dt * o[i](time,o_old[i]); for i from 1 to n do k[2,i]:=dt*o[i](time + dt/(2), o_old[i]+ (k[1, i])/(2)); end do: for i from 1 to n do k[3,i]:=dt*o[i](time+ o_old[i]+ ( k[1, i])/(2)); end do: for i from 1 to n do k[4,i]:=dt*o[i](time+dt,o_old[i]+k[1,i]); end do: for i from 1 to n do 0 _new[i]:=o_old[i]+(k[1,i]+2 *k[2,i]+2*k[3,i]+k[4,i]); end do: time := time + dt; end do; "

Error, illegal use of an object as a name

"local o_old; o_old := Vector(n); #` Loop to initialize o[i]old values` for i from 1 to n do o_old[i] := initial_conditions[i]; printf("Initial condition for u_%d: u_%d(0) = %f\\n", i, i, initial_conditions[i]); end do: time := t0; while time < tf do for i from 1 to n do k[1, i] := dt * o[i](time,o_old[i]); end do: for i from 1 to n do k[2,i]:=dt*o[i](time + dt/2, o_old[i]+ (k[1, i])/2); end do: for i from 1 to n do k[3,i]:=dt*o[i](time+ o_old[i]+ ( k[1, i])/2); end do: for i from 1 to n do k[4,i]:=dt*o[i](time+dt,o_old[i]+k[1,i]); end do: for i from 1 to n do 0 _new[i]:=o_old[i]+(k[1,i]+2 *k[2,i]+2*k[3,i]+k[4,i]); end do: time := time + dt; end do; "

  
 

NULL

Download dq_code.mw
Please Wait...