Advanced Search

janhardo

940 Reputation

13 Badges

12 years, 27 days
B. Ed math
Contact janhardo

MaplePrimes Activity


These are replies submitted by janhardo

corrected rod boundry conditions  ...

janhardo 940
7 hours ago
0 0

corrected rod boundry conditions 


 

> 

 

 

 
> 

restart:

interface(showassumed = 0):

with(plots):

#####################################################################
# EDUCATIONAL MAPLE WORKSHEET
# TRANSIENT AXISYMMETRIC HEAT CONDUCTION IN A FINITE CYLINDRICAL ROD
#
# Purpose
# -------
# This worksheet explains, derives, evaluates, and checks a double-series
# solution for a finite cylindrical rod subject to a radially varying,
# time-dependent heat flux q(r,t) at z = 0.
#
# The solution is axisymmetric:
#
#                       T = T(r,z,t)
#
# There is no dependence on the angular coordinate theta.
#
# IMPORTANT MAPLE NOTE
# --------------------
# The name gamma is protected in Maple.  Therefore this worksheet uses
# the name axialWaveNumber instead of gamma.
#
# Also note:
#     local declarations are used only inside proc ... end proc.
# All definitions outside procedures are global worksheet assignments.
#####################################################################


#####################################################################
# 0. PHYSICAL MODEL AND DOMAIN
#
# Domain:
#
#              0 <= r <= a
#              0 <= z <= L
#              t >= 0
#
# Variables:
#
# r : radial coordinate
# z : axial coordinate
# t : time
#
# Material constants:
#
# K     : thermal conductivity
# rho   : density
# cp    : specific heat capacity
# alpha : thermal diffusivity
#
# Relation:
#
#                   alpha = K/(rho*cp)
#
# Governing equation:
#
# rho*cp*T_t = K*(T_rr + (1/r)*T_r + T_zz)
#
# or equivalently:
#
# T_t = alpha*(T_rr + (1/r)*T_r + T_zz)
#
# Boundary conditions:
#
# BC1  -K*T_z(r,0,t) = q(r,t)      applied heat flux
# BC2      T_z(r,L,t) = 0           insulated end
# BC3        T(a,z,t) = 0           prescribed wall temperature
# BC4      T_r(0,z,t) = 0           axisymmetry
#
# Initial condition:
#
# IC         T(r,z,0) = 0
#
# CORNER COMPATIBILITY
# --------------------
# BC1 and BC3 meet at (r,z)=(a,0).
#
# Because BC3 gives T(a,z,t)=0 for every z, its axial derivative at the
# corner is also zero:
#
#                       T_z(a,0,t)=0.
#
# Then BC1 requires:
#
#                      q(a,t)=0.
#
# Hence q(a,t)=0 is necessary for classical compatibility.
#####################################################################

#####################################################################
# 1. SYMBOLIC PARAMETERS
#####################################################################

assume(K > 0):
assume(rho > 0):
assume(cp > 0):
assume(alpha > 0):
assume(a > 0):
assume(L > 0):

alpha_relation := alpha = K/(rho*cp):


#####################################################################
# 2. PDE, BOUNDARY CONDITIONS, AND INITIAL CONDITION
#####################################################################

PDE_cylinder :=
    diff(T(r,z,t),t)
    =
    alpha*(
        diff(T(r,z,t),r$2)
        +
        diff(T(r,z,t),r)/r
        +
        diff(T(r,z,t),z$2)
    ):

BC1_flux :=
    -K*D[2](T)(r,0,t)
    =
    q(r,t):

BC2_top :=
    D[2](T)(r,L,t)
    =
    0:

BC3_wall :=
    T(a,z,t)
    =
    0:

BC4_axis :=
    D[1](T)(0,z,t)
    =
    0:

IC :=
    T(r,z,0)
    =
    0:

corner_compatibility :=
    q(a,t)
    =
    0:


#####################################################################
# 3. RADIAL STURM-LIOUVILLE SYSTEM
#
# WHY BESSEL FUNCTIONS APPEAR
# ----------------------------
# Separation of variables in cylindrical coordinates produces
#
#        R'' + (1/r)R' + beta^2 R = 0.
#
# The solution regular at r=0 is J_0(beta*r).  The wall condition
# R(a)=0 selects discrete values beta_m=mu_m/a, where mu_m is the
# m-th positive zero of J_0.
#
# The radial eigenproblem is
#
# R'' + (1/r)R' + beta^2 R = 0,
# R'(0) = 0,
# R(a)  = 0.
#
# Hence:
#
# R_m(r) = J_0(mu_m*r/a),
#
# where mu_m is the m-th positive zero of J_0.
#
# beta_m = mu_m/a.
#
# Orthogonality:
#
# int(r*R_m(r)*R_j(r),r=0..a) = 0, m<>j,
#
# and
#
# int(r*R_m(r)^2,r=0..a)
#     = a^2*J_1(mu_m)^2/2.
#####################################################################

mu :=
    m ->
    BesselJZeros(0,m):

beta :=
    m ->
    mu(m)/a:

Rmode :=
    (m,rvalue) ->
    BesselJ(0,mu(m)*rvalue/a):

Rnorm :=
    m ->
    a^2*BesselJ(1,mu(m))^2/2:


#####################################################################
# 4. RADIAL EXPANSION OF THE APPLIED HEAT FLUX
#
# The flux profile must be expanded in the same radial eigenfunctions.
# The weight r in the integral is essential: it comes from the
# cylindrical area element and the Sturm-Liouville inner product.
#
# q(r,t) = Sum(q_m(t)*R_m(r),m=1..infinity),
#
# with
#
# q_m(t) =
#   2/(a^2*J_1(mu_m)^2)
#   * int(r*q(r,t)*J_0(mu_m*r/a),r=0..a).
#####################################################################

qmode :=
    (m,tvalue) ->
    2
    /
    (
        a^2*BesselJ(1,mu(m))^2
    )
    *
    Int(
        r*q(r,tvalue)*BesselJ(0,mu(m)*r/a),
        r=0..a
    ):


#####################################################################
# 5. AXIAL EIGENFUNCTIONS
#
# The axial eigenfunctions are cosines because the homogenized axial
# problem has Neumann conditions at z=0 and z=L.  The n=0 mode is the
# constant axial mode and must be retained.
#
# The homogeneous axial conditions are Neumann conditions:
#
# Z_n'(0)=0, Z_n'(L)=0.
#
# Thus:
#
# Z_0(z)=1,
# Z_n(z)=cos(n*Pi*z/L), n>=1.
#####################################################################

Zmode :=
    (n,zvalue) ->
    piecewise(
        n=0,
        1,
        cos(n*Pi*zvalue/L)
    ):

axialWaveNumber :=
    n ->
    n*Pi/L:

axial_factor :=
    n ->
    piecewise(
        n=0,
        1,
        2
    ):


#####################################################################
# 6. MODAL EQUATIONS
#
# Projection converts the PDE into independent first-order ODEs for
# the modal amplitudes A[m,n](t).  Each mode decays at a rate determined
# by radial diffusion plus axial diffusion.
#
# Write
#
# T(r,z,t)
# =
# Sum(
#   Sum(A[m,n](t)*R_m(r)*Z_n(z),n=0..infinity),
#   m=1..infinity
# ).
#
# Projection gives
#
# A[m,n]'(t)
# + alpha*(beta_m^2 + gamma_n^2)*A[m,n](t)
#
# =
# axial_factor(n)*alpha*q_m(t)/(K*L),
#
# with A[m,n](0)=0.
#
# Therefore:
#
# A[m,n](t)
# =
# axial_factor(n)*alpha/(K*L)
# *
# int(
#   exp(-alpha*(beta_m^2+gamma_n^2)*(t-s))*q_m(s),
#   s=0..t
# ).
#####################################################################

decay :=
    (m,n) ->
    alpha*(
        beta(m)^2
        +
        axialWaveNumber(n)^2
    ):

Amn :=
    (m,n,tvalue) ->
    axial_factor(n)*alpha/(K*L)
    *
    Int(
        exp(
            -decay(m,n)*(tvalue-s)
        )
        *
        qmode(m,s),
        s=0..tvalue
    ):


#####################################################################
# 7. FORMAL DOUBLE-SERIES SOLUTION
#
# This is a formal symbolic representation.  Maple may display inert
# Int and Sum objects.  Numerical work below uses finite sums and
# evaluates the radial projections numerically.
#####################################################################

T_exact :=
    Sum(
        Sum(
            Amn(m,n,t)
            *
            Rmode(m,r)
            *
            Zmode(n,z),
            n=0..infinity
        ),
        m=1..infinity
    ):

printf("\nFormal double-series solution:\n"):

T(r,z,t) = T_exact;


#####################################################################
# 8. A COMPATIBLE EXAMPLE FLUX
#
# The factor 1-(r/a)^2 forces q(a,t)=0.  The factor sin(omega*t)
# forces q(r,0)=0, making the flux compatible with the zero initial
# temperature at t=0.
#
# Choose
#
# q(r,t) = q0*(1-(r/a)^2)*sin(omega*t).
#
# This satisfies:
#
# q(a,t)=0 for every t,
#
# and also q(r,0)=0, so the flux is compatible with the zero initial
# temperature at t=0.
#####################################################################

assume(q0 > 0):
assume(omega > 0):

q_example :=
    (rvalue,tvalue) ->
    q0*(
        1-(rvalue/a)^2
    )
    *
    sin(omega*tvalue):

corner_check :=
    simplify(
        q_example(a,t)
    ):

printf(
    "\nCorner compatibility check q(a,t): %a\n",
    corner_check
):

initial_flux_check :=
    simplify(
        q_example(r,0)
    ):

printf(
    "Initial flux compatibility check q(r,0): %a\n",
    initial_flux_check
):


#####################################################################
# 9. RADIAL COEFFICIENTS FOR THE EXAMPLE FLUX
#
# Since q(r,t)=f(r)*sin(omega*t), write
#
# q_m(t)=Q_m*sin(omega*t).
#####################################################################

Qm :=
    m ->
    simplify(
        2*q0
        /
        (
            a^2*BesselJ(1,mu(m))^2
        )
        *
        Int(
            r
            *
            (
                1-(r/a)^2
            )
            *
            BesselJ(0,mu(m)*r/a),
            r=0..a
        )
    ):

qmode_example :=
    (m,tvalue) ->
    Qm(m)*sin(omega*tvalue):


#####################################################################
# 10. CLOSED TIME CONVOLUTION FOR SINUSOIDAL HEATING
#
# For lambda_mn = alpha*(beta_m^2+gamma_n^2),
#
# int(
#   exp(-lambda_mn*(t-s))*sin(omega*s),
#   s=0..t
# )
#
# =
# (
#   lambda_mn*sin(omega*t)
#   -omega*cos(omega*t)
#   +omega*exp(-lambda_mn*t)
# )
# /
# (
#   lambda_mn^2+omega^2
# ).
#####################################################################

TimeConvolution :=
    (m,n,tvalue) ->
    (
        decay(m,n)*sin(omega*tvalue)
        -
        omega*cos(omega*tvalue)
        +
        omega*exp(-decay(m,n)*tvalue)
    )
    /
    (
        decay(m,n)^2
        +
        omega^2
    ):

Amn_example :=
    (m,n,tvalue) ->
    axial_factor(n)*alpha/(K*L)
    *
    Qm(m)
    *
    TimeConvolution(m,n,tvalue):

T_example_exact :=
    Sum(
        Sum(
            Amn_example(m,n,t)
            *
            Rmode(m,r)
            *
            Zmode(n,z),
            n=0..infinity
        ),
        m=1..infinity
    ):

printf(
    "\nSeries solution for q(r,t)=q0*(1-(r/a)^2)*sin(omega*t):\n"
):

T(r,z,t) = T_example_exact;


#####################################################################
# 11. NUMERICAL PARAMETER VALUES
#
# Replace these demonstration values by the physical values of the rod.
# Units must be consistent.  In SI units:
# K [W/(m K)], rho [kg/m^3], cp [J/(kg K)], a,L [m], q0 [W/m^2].
#
# K     : thermal conductivity
# rho   : density
# cp    : specific heat
# alpha : thermal diffusivity K/(rho*cp)
#####################################################################

Kval :=
    10.0:

rhoval :=
    7800.0:

cpval :=
    500.0:

alphaval :=
    Kval/(rhoval*cpval):

aval :=
    0.5:

Lval :=
    2.0:

q0val :=
    1000.0:

omegaval :=
    1.0:

Mterms :=
    15:

Nterms :=
    25:


#####################################################################
# 12. NUMERICAL RADIAL EIGENVALUES AND EIGENFUNCTIONS
#
# Bessel zeros are cached implicitly by Maple when repeatedly evaluated.
# The radial eigenfunctions automatically satisfy T(a,z,t)=0.
#####################################################################

mu_num :=
    m ->
    evalf(
        BesselJZeros(0,m)
    ):

beta_num :=
    m ->
    mu_num(m)/aval:

Rmode_num :=
    (m,rvalue) ->
    BesselJ(
        0,
        mu_num(m)*rvalue/aval
    ):

gamma_num :=
    n ->
    n*Pi/Lval:

axial_factor_num :=
    n ->
    piecewise(
        n=0,
        1.0,
        2.0
    ):

decay_num :=
    (m,n) ->
    alphaval*(
        beta_num(m)^2
        +
        gamma_num(n)^2
    ):


#####################################################################
# 13. NUMERICAL RADIAL FLUX COEFFICIENTS
#
# Qm_num is a Maple procedure because it needs a local variable and
# option remember.  This avoids recalculating the same quadrature for
# every plot point.
#
# Hardware quadrature is used once for each radial mode.
#####################################################################

Qm_num :=
proc(m)

    option remember;

    local mum;

    mum :=
        mu_num(m);

    return evalf(
        2*q0val
        /
        (
            aval^2*BesselJ(1,mum)^2
        )
        *
        evalf(
            Int(
                x
                *
                (
                    1-(x/aval)^2
                )
                *
                BesselJ(
                    0,
                    mum*x/aval
                ),
                x=0..aval,
                method=_Gquad
            )
        )
    );

end proc:


#####################################################################
# 14. NUMERICAL MODAL AMPLITUDES
#####################################################################

TimeConvolution_num :=
    (m,n,tvalue) ->
    (
        decay_num(m,n)*sin(omegaval*tvalue)
        -
        omegaval*cos(omegaval*tvalue)
        +
        omegaval*exp(-decay_num(m,n)*tvalue)
    )
    /
    (
        decay_num(m,n)^2
        +
        omegaval^2
    ):

Amn_num :=
    (m,n,tvalue) ->
    axial_factor_num(n)
    *
    alphaval
    /
    (
        Kval*Lval
    )
    *
    Qm_num(m)
    *
    TimeConvolution_num(m,n,tvalue):


#####################################################################
# 15. TRUNCATED DOUBLE SERIES
#
# The exact solution uses infinitely many radial and axial modes.
# Numerically, we truncate at m=1..M and n=0..N.
# Increase M and N until the result changes negligibly.
#####################################################################

Tnum :=
proc(rvalue,zvalue,tvalue,M,N)

    local m;
    local n;
    local total;

    if not type(rvalue,numeric) then
        return 'procname'(rvalue,zvalue,tvalue,M,N);
    end if;

    if not type(zvalue,numeric) then
        return 'procname'(rvalue,zvalue,tvalue,M,N);
    end if;

    if not type(tvalue,numeric) then
        return 'procname'(rvalue,zvalue,tvalue,M,N);
    end if;

    total :=
        add(
            add(
                Amn_num(m,n,tvalue)
                *
                Rmode_num(m,rvalue)
                *
                piecewise(
                    n=0,
                    1.0,
                    cos(n*Pi*zvalue/Lval)
                ),
                n=0..N
            ),
            m=1..M
        );

    return evalf(total);

end proc:


#####################################################################
# 16. NUMERICAL CHECKS
#
# A small nonzero wall value, for example 1e-14, is numerical roundoff
# and should be interpreted as zero.
#####################################################################

printf(
    "\nNumerical q(a,t) check at t=1: %g\n",
    q0val*(1-(aval/aval)^2)*sin(omegaval)
):

printf(
    "Initial temperature check T(r,z,0): %g\n",
    Tnum(0.2,0.7,0.0,Mterms,Nterms)
):

printf(
    "Outer-wall check T(a,z,t): %g\n",
    Tnum(aval,0.7,1.0,Mterms,Nterms)
):

printf(
    "Sample temperature T(0.2,0.7,1): %g\n",
    Tnum(0.2,0.7,1.0,Mterms,Nterms)
):


#####################################################################
# 17. TEMPERATURE ON THE r-z CROSS-SECTION AT A FIXED TIME
#####################################################################

tplot :=
    2.0:

CrossSectionTemperature :=
plot3d(
    (rvalue,zvalue) ->
        Tnum(
            rvalue,
            zvalue,
            tplot,
            Mterms,
            Nterms
        ),

    0..aval,
    0..Lval,

    labels =
    [
        "radius r",
        "axial coordinate z",
        "temperature T(r,z,t)"
    ],

    title =
    cat(
        "Axisymmetric temperature field at t = ",
        tplot
    ),

    axes =
    boxed,

    grid =
    [35,45],

    orientation =
    [55,65]
):

CrossSectionTemperature;


#####################################################################
# 18. RADIAL TEMPERATURE PROFILES AT z=0
#####################################################################

RadialProfiles :=
plot(
    [
        rvalue ->
            Tnum(
                rvalue,
                0.0,
                0.5,
                Mterms,
                Nterms
            ),

        rvalue ->
            Tnum(
                rvalue,
                0.0,
                1.0,
                Mterms,
                Nterms
            ),

        rvalue ->
            Tnum(
                rvalue,
                0.0,
                2.0,
                Mterms,
                Nterms
            )
    ],

    0..aval,

    labels =
    [
        "radius r",
        "temperature T(r,0,t)"
    ],

    legend =
    [
        "t = 0.5",
        "t = 1.0",
        "t = 2.0"
    ],

    title =
    "Radial temperature profiles on the heated end",

    thickness =
    2,

    adaptive =
    false,

    numpoints =
    180
):

RadialProfiles;




#####################################################################
# 19. EDUCATIONAL CONVERGENCE STUDY
#
# A truncated series should be checked by increasing M and N.
# The following procedure compares several truncation levels at one
# selected point (rtest,ztest,ttest).
#####################################################################

ConvergenceStudy :=
proc(rtest,ztest,ttest)

    local truncations;
    local j;
    local Mj;
    local Nj;
    local valuej;

    truncations :=
        [
            [5,8],
            [8,12],
            [12,18],
            [15,25],
            [20,35]
        ];

    printf(
        "\nConvergence study at r=%g, z=%g, t=%g\n",
        rtest,
        ztest,
        ttest
    );

    printf(
        "      M       N              T_MN\n"
    );

    for j from 1 to nops(truncations) do

        Mj := truncations[j][1];
        Nj := truncations[j][2];

        valuej :=
            Tnum(
                rtest,
                ztest,
                ttest,
                Mj,
                Nj
            );

        printf(
            "%7d %7d %20.12g\n",
            Mj,
            Nj,
            valuej
        );

    end do;

    return NULL;

end proc:

ConvergenceStudy(
    0.2,
    0.7,
    1.0
):


#####################################################################
# 20. EDUCATIONAL SUMMARY
#
# 1. Because q=q(r,t), radial dependence cannot be discarded.
#
# 2. The wall condition T(a,z,t)=0 gives Bessel J_0 radial modes.
#
# 3. The insulated axial end and homogenized lower-end condition give
#    cosine modes in z.
#
# 4. The compatibility condition q(a,t)=0 is required at the corner
#    where BC1 and BC3 meet.
#
# 5. The exact solution is a double infinite series.
#
# 6. Tnum evaluates a truncated version of that series.
#
# 7. Mterms and Nterms should be selected through a convergence study.
#####################################################################

printf(
    "\nWorksheet completed successfully.\n"
):

#####################################################################
# END OF THE EDUCATIONAL MAPLE WORKSHEET
#####################################################################


Formal double-series solution:

  

T(r, z, t) = Sum(Sum(piecewise(n = 0, 1, 2)*alpha*(Int(2*exp(-alpha*(BesselJZeros(0, m)^2/a^2+n^2*Pi^2/L^2)*(t-s))*(Int(r*q(r, s)*BesselJ(0, BesselJZeros(0, m)*r/a), r = 0 .. a))/(a^2*BesselJ(1, BesselJZeros(0, m))^2), s = 0 .. t))*BesselJ(0, BesselJZeros(0, m)*r/a)*piecewise(n = 0, 1, cos(n*Pi*z/L))/(K*L), n = 0 .. infinity), m = 1 .. infinity)

  


Corner compatibility check q(a,t): 0
Initial flux compatibility check q(r,0): 0

Series solution for q(r,t)=q0*(1-(r/a)^2)*sin(omega*t):

  

T(r, z, t) = Sum(Sum(2*piecewise(n = 0, 1, 2)*alpha*q0*(Int(r*(a^2-r^2)*BesselJ(0, BesselJZeros(0, m)*r/a), r = 0 .. a))*(alpha*(BesselJZeros(0, m)^2/a^2+n^2*Pi^2/L^2)*sin(omega*t)-omega*cos(omega*t)+omega*exp(-alpha*(BesselJZeros(0, m)^2/a^2+n^2*Pi^2/L^2)*t))*BesselJ(0, BesselJZeros(0, m)*r/a)*piecewise(n = 0, 1, cos(n*Pi*z/L))/(K*L*a^4*BesselJ(1, BesselJZeros(0, m))^2*(alpha^2*(BesselJZeros(0, m)^2/a^2+n^2*Pi^2/L^2)^2+omega^2)), n = 0 .. infinity), m = 1 .. infinity)

  


Numerical q(a,t) check at t=1: 0
Initial temperature check T(r,z,0): 0
Outer-wall check T(a,z,t): -4.20561e-15
Sample temperature T(0.2,0.7,1): 2.20970e-05

  

 

 


Convergence study at r=0.2, z=0.7, t=1
      M       N              T_MN
      5       8      7.449060185e-06
      8      12      8.754690157e-05
     12      18      9.439198208e-05
     15      25      2.209696623e-05
     20      35      9.187070509e-05

Worksheet completed successfully.

  
> 

 

> 

 


 

Download pde_heat_1D_oplossing_van_mprimes_misschien_correctie_17-7-2026.mw

 

boundry conditions for this rod setup&nb...

janhardo 940
8 hours ago
0 0

boundry conditions for this rod setup 

@Rouben Rostamian  I'm not enti...

janhardo 940
9 hours ago
0 0

@Rouben Rostamian  
Thanks for pointing this out
I'm not entirely sure what you mean exactly, but got a idea of it 

@C_R I tried something different, b...

janhardo 940
11 hours ago
0 0

@C_R 
I tried something different, but not yet familiar with techniques used here 

> 

 

pde heat 1 D

 
> 

restart:

interface(showassumed = 0):

with(plots):

#####################################################################
# TRANSIENT HEAT CONDUCTION IN AN INSULATED CYLINDRICAL ROD
#
# Original cylindrical problem:
#
# T_t = k*(T_rr + (1/r)*T_r + T_zz)
#
# Boundary conditions:
#
# -k*T_z(r,0,t) = q(t)       prescribed heat flux at z = 0
#     T_z(r,L,t) = 0         insulated upper end
#     T_r(R,z,t) = 0         insulated outer wall
#     T_r(0,z,t) = 0         symmetry at the cylinder axis
#
# Initial condition:
#
# T(r,z,0) = 0
#
# Since all data are independent of r and the outer wall is insulated,
# the temperature is radially uniform:
#
#                   T(r,z,t) = V(z,t)
#
# Therefore:
#
# T_r = 0,  T_rr = 0
#
# and the cylindrical problem reduces to:
#
# V_t = k*V_zz
#
# -k*V_z(0,t) = q(t)
#     V_z(L,t) = 0
#     V(z,0)   = 0
#####################################################################


#####################################################################
# 1. SYMBOLIC PARAMETERS
#####################################################################

assume(k > 0):
assume(L > 0):
assume(R > 0):


#####################################################################
# 2. ORIGINAL CYLINDRICAL PDE
#####################################################################

PDE_cylinder :=
    diff(T(r,z,t),t)
    =
    k*(
        diff(T(r,z,t),r$2)
        +
        diff(T(r,z,t),r)/r
        +
        diff(T(r,z,t),z$2)
    ):

BC_cylinder_bottom :=
    -k*D[2](T)(r,0,t)
    =
    q(t):

BC_cylinder_top :=
    D[2](T)(r,L,t)
    =
    0:

BC_cylinder_wall :=
    D[1](T)(R,z,t)
    =
    0:

BC_cylinder_axis :=
    D[1](T)(0,z,t)
    =
    0:

IC_cylinder :=
    T(r,z,0)
    =
    0:


#####################################################################
# 3. REDUCED ONE-DIMENSIONAL PROBLEM
#
# T(r,z,t) = V(z,t)
#####################################################################

PDE_1D :=
    diff(V(z,t),t)
    =
    k*diff(V(z,t),z$2):

BC_bottom :=
    -k*D[1](V)(0,t)
    =
    q(t):

BC_top :=
    D[1](V)(L,t)
    =
    0:

IC_1D :=
    V(z,0)
    =
    0:

problem_1D :=
    {
        PDE_1D,
        BC_bottom,
        BC_top,
        IC_1D
    }:


#####################################################################
# 4. STURM-LIOUVILLE EIGENVALUE PROBLEM
#
# X''(z) + lambda*X(z) = 0
#
# X'(0) = 0
# X'(L) = 0
#
# The eigenfunctions are:
#
# phi_0(z) = 1
#
# phi_n(z) = cos(n*Pi*z/L),   n = 1,2,3,...
#
# The eigenvalues are:
#
# lambda_0 = 0
#
# lambda_n = (n*Pi/L)^2
#####################################################################

phi0 :=
    z -> 1:

phi :=
    (n,z) ->
    cos(n*Pi*z/L):

lambda :=
    n ->
    (n*Pi/L)^2:


#####################################################################
# 5. VERIFY THE STURM-LIOUVILLE EQUATION
#
# phi_n'' + lambda_n*phi_n = 0
#####################################################################

SL_equation_check :=
    simplify(
        diff(phi(n,z),z$2)
        +
        lambda(n)*phi(n,z)
    ):

printf(
    "Sturm-Liouville differential equation check: %a\n",
    SL_equation_check
):


#####################################################################
# 6. VERIFY THE HOMOGENEOUS NEUMANN CONDITIONS
#####################################################################

SL_left_BC_check :=
    simplify(
        subs(
            z=0,
            diff(phi(n,z),z)
        )
    ):

SL_right_BC_check :=
    simplify(
        subs(
            z=L,
            diff(phi(n,z),z)
        )
    ):

printf(
    "Eigenfunction derivative at z=0: %a\n",
    SL_left_BC_check
):

printf(
    "Eigenfunction derivative at z=L: %a\n",
    SL_right_BC_check
):


#####################################################################
# 7. NORMS OF THE EIGENFUNCTIONS
#
# Integral(phi_0^2,z=0..L) = L
#
# Integral(phi_n^2,z=0..L) = L/2, n >= 1
#####################################################################

Norm0 :=
    int(
        phi0(z)^2,
        z=0..L
    ):

NormN_formula :=
    L/2:

printf(
    "Norm squared of the zero mode: %a\n",
    Norm0
):

printf(
    "Norm squared of every mode n>=1: %a\n",
    NormN_formula
):


#####################################################################
# 8. FOURIER-COSINE REPRESENTATION
#
# V(z,t) =
#
# a0(t)
# +
# Sum(
#     a_n(t)*cos(n*Pi*z/L),
#     n=1..infinity
# )
#####################################################################

V_Fourier :=
    a0(t)
    +
    Sum(
        a(n,t)*phi(n,z),
        n=1..infinity
    ):


#####################################################################
# 9. ZERO-MODE EQUATION
#
# Integrate the PDE over 0 < z < L:
#
# Integral(V_t dz)
# =
# k*[V_z(L,t)-V_z(0,t)]
#
# Because:
#
# V_z(L,t) = 0
#
# V_z(0,t) = -q(t)/k
#
# we obtain:
#
# d/dt Integral(V dz) = q(t)
#
# Since:
#
# a0(t) = (1/L)*Integral(V dz)
#
# the zero-mode equation is:
#
# a0'(t) = q(t)/L
#
# a0(0) = 0
#####################################################################

ODE_zero_mode :=
    diff(a0(t),t)
    =
    q(t)/L:

IC_zero_mode :=
    a0(0)
    =
    0:

a0_solution :=
    Int(
        q(s),
        s=0..t
    )/L:

zero_mode_check :=
    simplify(
        diff(a0_solution,t)
        -
        q(t)/L
    ):

printf(
    "Zero-mode ODE check: %a\n",
    zero_mode_check
):


#####################################################################
# 10. HIGHER-MODE EQUATIONS
#
# Multiply the PDE by phi_n(z) and integrate from 0 to L.
#
# After two integrations by parts:
#
# a_n'(t) + k*lambda_n*a_n(t) = 2*q(t)/L
#
# with:
#
# a_n(0) = 0
#####################################################################

ODE_mode_n :=
    diff(a(n,t),t)
    +
    k*lambda(n)*a(n,t)
    =
    2*q(t)/L:

IC_mode_n :=
    a(n,0)
    =
    0:


#####################################################################
# 11. SOLUTION OF EACH HIGHER-MODE ODE
#
# Using the integrating factor:
#
# exp(k*lambda_n*t)
#
# we obtain:
#
# a_n(t) =
#
# (2/L)*Integral(
#     exp(-k*lambda_n*(t-s))*q(s),
#     s=0..t
# )
#####################################################################

an_solution :=
    (2/L)*
    Int(
        exp(
            -k*lambda(n)*(t-s)
        )*
        q(s),
        s=0..t
    ):

mode_n_check :=
    simplify(
        diff(an_solution,t)
        +
        k*lambda(n)*an_solution
        -
        2*q(t)/L
    ):

printf(
    "Higher-mode ODE check: %a\n",
    mode_n_check
):


#####################################################################
# 12. GENERAL ANALYTICAL SOLUTION
#
# V(z,t) =
#
# (1/L)*Integral(q(s),s=0..t)
#
# +
#
# (2/L)*Sum(
#     cos(n*Pi*z/L)
#     *
#     Integral(
#         exp(-k*(n*Pi/L)^2*(t-s))*q(s),
#         s=0..t
#     ),
#     n=1..infinity
# )
#####################################################################

V_exact :=
    Int(
        q(s),
        s=0..t
    )/L
    +
    (2/L)*
    Sum(
        cos(n*Pi*z/L)
        *
        Int(
            exp(
                -k*(n*Pi/L)^2*(t-s)
            )*
            q(s),
            s=0..t
        ),
        n=1..infinity
    ):

printf(
    "\nGeneral analytical solution:\n"
):

V(z,t) = V_exact;


#####################################################################
# 13. RETURN TO THE ORIGINAL CYLINDRICAL ROD
#
# Because there is no radial temperature variation:
#
# T(r,z,t) = V(z,t)
#####################################################################

T_exact :=
    V_exact:

radial_derivative_check :=
    simplify(
        diff(T_exact,r)
    ):

printf(
    "Radial derivative T_r = %a\n",
    radial_derivative_check
):


#####################################################################
# 14. INITIAL-CONDITION CHECK
#####################################################################

initial_condition_check :=
    simplify(
        subs(
            t=0,
            V_exact
        )
    ):

printf(
    "Initial-condition check V(z,0) = %a\n",
    initial_condition_check
):


#####################################################################
# 15. SPECIFIC HEAT FLUX
#
# Choose:
#
# q(t) = q0*sin(omega*t)
#
# This is convenient because:
#
# q(0) = 0
#
# so the initial and boundary data are compatible at t=0.
#####################################################################

assume(q0 > 0):
assume(omega > 0):

q_sinusoidal :=
    tau ->
    q0*sin(omega*tau):


#####################################################################
# 16. ZERO MODE FOR THE SINUSOIDAL FLUX
#
# a0(t) =
#
# q0*(1-cos(omega*t))/(L*omega)
#####################################################################

a0_sinusoidal :=
    simplify(
        int(
            q_sinusoidal(s),
            s=0..t
        )/L
    ):

printf(
    "\nZero mode for q(t)=q0*sin(omega*t):\n"
):

a0(t) = a0_sinusoidal;


#####################################################################
# 17. CLOSED FORM OF EACH HIGHER MODE
#
# Define:
#
# alpha_n = k*(n*Pi/L)^2
#
# Then:
#
# Integral(
#   exp(-alpha_n*(t-s))*q0*sin(omega*s),
#   s=0..t
# )
#
# =
#
# q0*
# (
#   alpha_n*sin(omega*t)
#   -omega*cos(omega*t)
#   +omega*exp(-alpha_n*t)
# )
# /
# (
#   alpha_n^2 + omega^2
# )
#####################################################################

alpha_symbolic :=
    n ->
    k*(n*Pi/L)^2:

convolution_sinusoidal :=
    (n,t) ->
    q0*
    (
        alpha_symbolic(n)*sin(omega*t)
        -
        omega*cos(omega*t)
        +
        omega*exp(-alpha_symbolic(n)*t)
    )
    /
    (
        alpha_symbolic(n)^2
        +
        omega^2
    ):


#####################################################################
# 18. EXACT SERIES FOR THE SINUSOIDAL FLUX
#####################################################################

V_sinusoidal :=
    q0*
    (
        1-cos(omega*t)
    )
    /
    (
        L*omega
    )
    +
    (2/L)*
    Sum(
        cos(n*Pi*z/L)
        *
        convolution_sinusoidal(n,t),
        n=1..infinity
    ):

printf(
    "\nAnalytical Fourier-series solution for the sinusoidal flux:\n"
):

V(z,t) = V_sinusoidal;


#####################################################################
# 19. NUMERICAL VALUES
#
# Replace these values by the physical values of the experiment.
#####################################################################

kval :=
    0.1:

Lval :=
    2.0:

Rval :=
    0.5:

q0val :=
    10.0:

omegaval :=
    1.0:

Nterms :=
    40:


#####################################################################
# 20. NUMERICAL DECAY CONSTANT
#
# alpha_n = k*(n*Pi/L)^2
#####################################################################

alpha_num :=
    n ->
    kval*(n*Pi/Lval)^2:


#####################################################################
# 21. NUMERICAL ZERO MODE
#####################################################################

MeanTerm_num :=
    tvalue ->
    q0val*
    (
        1-cos(omegaval*tvalue)
    )
    /
    (
        Lval*omegaval
    ):


#####################################################################
# 22. NUMERICAL MODAL CONVOLUTION
#
# The time integral has already been evaluated analytically.
#
# This is essential for fast plotting.
#####################################################################

ModalConvolution_num :=
    (n,tvalue) ->
    q0val*
    (
        alpha_num(n)*sin(omegaval*tvalue)
        -
        omegaval*cos(omegaval*tvalue)
        +
        omegaval*exp(-alpha_num(n)*tvalue)
    )
    /
    (
        alpha_num(n)^2
        +
        omegaval^2
    ):


#####################################################################
# 23. TRUNCATED FOURIER SERIES
#
# N controls the number of Fourier modes.
#
# Larger N gives a more accurate approximation but requires
# more computation.
#####################################################################

Vnum :=
proc(zvalue,tvalue,N)

    local n;
    local MeanPart;
    local SeriesPart;

    if not type(zvalue,numeric) then
        return 'procname'(zvalue,tvalue,N);
    end if;

    if not type(tvalue,numeric) then
        return 'procname'(zvalue,tvalue,N);
    end if;

    MeanPart :=
        MeanTerm_num(tvalue);

    SeriesPart :=
        (2/Lval)*
        add(
            cos(n*Pi*zvalue/Lval)
            *
            ModalConvolution_num(n,tvalue),
            n=1..N
        );

    return evalf(
        MeanPart
        +
        SeriesPart
    );

end proc:


#####################################################################
# 24. TEST THE NUMERICAL FUNCTION BEFORE PLOTTING
#####################################################################

TestTemperature :=
    Vnum(
        0.5,
        3.0,
        Nterms
    ):

printf(
    "\nTemperature at z=0.5 and t=3:\n"
):

TestTemperature;


#####################################################################
# 25. CHECK THE INITIAL CONDITION NUMERICALLY
#####################################################################

InitialValueAtTestPoint :=
    Vnum(
        0.5,
        0.0,
        Nterms
    ):

printf(
    "\nNumerical initial-condition check:\n"
):

InitialValueAtTestPoint;


#####################################################################
# 26. TEMPERATURE PROFILES AT DIFFERENT TIMES
#
# The plot is returned directly.
#
# There is no separate plots:-display call.
#####################################################################

TemperatureProfiles :=
plot(
    [
        zvalue -> Vnum(zvalue,0.25,Nterms),
        zvalue -> Vnum(zvalue,0.50,Nterms),
        zvalue -> Vnum(zvalue,1.00,Nterms),
        zvalue -> Vnum(zvalue,2.00,Nterms),
        zvalue -> Vnum(zvalue,3.00,Nterms)
    ],
    0..Lval,

    labels =
    [
        "axial coordinate z",
        "temperature V(z,t)"
    ],

    legend =
    [
        "t = 0.25",
        "t = 0.50",
        "t = 1.00",
        "t = 2.00",
        "t = 3.00"
    ],

    title =
    "Transient temperature profiles in the rod",

    thickness =
    2,

    adaptive =
    false,

    numpoints =
    180
):

TemperatureProfiles;


#####################################################################
# 27. THREE-DIMENSIONAL SPACE-TIME TEMPERATURE SURFACE
#
# For a faster first plot, use:
#
# Nterms = 20
# grid    = [25,25]
#
# For a smoother result, use:
#
# Nterms = 40
# grid    = [40,40]
#####################################################################

TemperatureSurface :=
plot3d(
    (zvalue,tvalue) ->
        Vnum(
            zvalue,
            tvalue,
            Nterms
        ),

    0..Lval,
    0..6,

    labels =
    [
        "axial coordinate z",
        "time t",
        "temperature V(z,t)"
    ],

    title =
    "Transient heat conduction in the insulated rod",

    axes =
    boxed,

    grid =
    [35,35],

    orientation =
    [55,65]
):

TemperatureSurface;


#####################################################################
# 28. OPTIONAL: TEMPERATURE ALONG THE ORIGINAL CYLINDER
#
# Since T(r,z,t)=V(z,t), the temperature is constant in r.
#
# This plot shows the same axial temperature profile repeated over
# the complete cylinder radius.
#####################################################################

CylinderTemperature :=
plot3d(
    (rvalue,zvalue) ->
        Vnum(
            zvalue,
            3.0,
            Nterms
        ),

    0..Rval,
    0..Lval,

    labels =
    [
        "radius r",
        "axial coordinate z",
        "temperature T(r,z,3)"
    ],

    title =
    "Temperature over the cylindrical cross-section at t=3",

    axes =
    boxed,

    grid =
    [25,40],

    orientation =
    [55,65]
):

CylinderTemperature;


#####################################################################
# END OF THE MAPLE WORKSHEET
#####################################################################

Sturm-Liouville differential equation check: 0
Eigenfunction derivative at z=0: 0
Eigenfunction derivative at z=L: -n*Pi/L*sin(n*Pi)
Norm squared of the zero mode: L
Norm squared of every mode n>=1: 1/2*L
Zero-mode ODE check: 0
Higher-mode ODE check: 0

General analytical solution:

  

V(z, t) = (Int(q(s), s = 0 .. t))/L+2*(Sum(cos(n*Pi*z/L)*(Int(exp(-k*n^2*Pi^2*(t-s)/L^2)*q(s), s = 0 .. t)), n = 1 .. infinity))/L

 

Radial derivative T_r = 0
Initial-condition check V(z,0) = 0

Zero mode for q(t)=q0*sin(omega*t):

  

a0(t) = -q0*(cos(omega*t)-1)/(omega*L)

 


Analytical Fourier-series solution for the sinusoidal flux:

  

V(z, t) = q0*(-cos(omega*t)+1)/(L*omega)+2*(Sum(cos(n*Pi*z/L)*q0*(k*n^2*Pi^2*sin(omega*t)/L^2-omega*cos(omega*t)+omega*exp(-k*n^2*Pi^2*t/L^2))/(k^2*n^4*Pi^4/L^4+omega^2), n = 1 .. infinity))/L

 


Temperature at z=0.5 and t=3:

  

17.38319022

 


Numerical initial-condition check:

  

0.

 

 

 

 
> 

 

Download pde_heat_1D_oplossing_van_mprimes_misschien_17-7-2026.mw

 restart:infolevel[pdsolve] := 3:as...

janhardo 940
16 hours ago
1 0

https://www.dropbox.com/scl/fi/cuexecaf3x4wd5vwi28gv/pde-warmet-cilindrische-staaf.png?rlkey=9gob68msfpsdbksuemcxtwtiz&st=t6d7q09s&dl=0

 

restart:

infolevel[pdsolve] := 3:

assume(k > 0, L > 0):

PDE1D :=
    diff(V(z,t),t)
    =
    k*diff(V(z,t),z$2):

BC1 :=
    -k*D[1](V)(0,t)
    =
    q(t):

BC2 :=
    D[1](V)(L,t)
    =
    0:

IC :=
    V(z,0)
    =
    0:

problem1D := {PDE1D, BC1, BC2, IC}:

sol1D := pdsolve(
    problem1D,
    V(z,t)

This seems the right conditions for this...

janhardo 940
18 hours ago
0 0

This seems the right conditions for this experiment ?

BC1 := -k*D[2](T)(r,0,t) = q(t);
BC2 := D[2](T)(r,L,t) = 0;
BC3 := D[1](T)(R,z,t) = 0;
BC4 := D[1](T)(0,z,t) = 0;
IC  := T(r,z,0) = 0;

 

Improved experiments heat_pde_proce...

janhardo 940
Yesterday at 5:14 PM
0 0

Improved experiments, but handling pdsolve command is something different 
 heat_pde_procedure_mprimes-echte_temperaturr_16-7-2026.mw                                                   

@Earl  Hello Earl , the display com...

janhardo 940
July 08 2026
0 0

@Earl 
Hello Earl , the display command ends with : (surpressing outcome),  use  ;  

3_regions_-_proceedure_-velocity_-temper...

janhardo 940
July 07 2026
0 0

3_regions_-_proceedure_-velocity_-temperatuur_profielen_7-7-2026.mw
a array plot ...?

Most close to get , so there is more involved ?
most_close_table_to_get_mprimes_7-7-2026.mw
Module code ( start)
3_regions_mprimes-verbeteringen_8-7-2026.mw

 

@dharr Ctr J, K  and Ctrl T gi...

janhardo 940
July 02 2026
0 0

@dharr 
Thanks, Ctrl J, K  and Ctrl T gives text  and Ctrl M / Ctrl R  seems to me the important key combinations for working in Maple. 
Ctrl .  gives a section , handy too.
Ctrl  Z and Ctrl  X , do and undo. 

Does it depend on the solution family of...

janhardo 940
June 22 2026
0 0

Does it depend on the solution family of the ODE that Maple apparently chooses?

Your expectation that the ODE output always uses set notation is probably not always possible.

@acer  "Just to be clear,...

janhardo 940
June 17 2026
0 0

@acer 

"Just to be clear, that P formula represents a product of terms (and not a sequence)."?

Just to be clear, that P formula represents a product of factors (and not a sequence)

Terms are used in a sum. 
Something went wrong with the message I wrote earlier, and a sentence was left out.
So it is a correction on your earlier sentence.

@acer "Just to be clear, that ...

janhardo 940
June 16 2026
0 0

@acer 
"Just to be clear, that P formula represents a product of terms (and not a sequence)."?

Terms are related to a sum.

@Alfred_F Earlier versions of ...

janhardo 940
June 04 2026
0 0

@Alfred_F 
Earlier versions of  Maple :  Maple's start-up message, 'Command the brilliance of 1,000 mathematicians'. :-)
 

So yes, from a mathematician's persp...

janhardo 940
June 03 2026
0 0

Thus, from a mathematician’s perspective, this framework is essentially a symbolic ansatz engine expressed in neural-network notation. This also explains why the method ultimately yields explicit analytical formulas, such as Equation (12), rather than a trained predictive model.

In other words, the neural-network architecture serves primarily as a structured language for generating candidate solution families, while the actual determination of the parameters is carried out analytically through algebraic constraints derived from the differential equation. The process therefore resembles symbolic computation and exact-solution techniques far more than modern machine-learning approaches based on optimization and training.

1 2 3 4 5 6 7 Last Page 1 of 88