Computational analysis of a class of singular nonlinear fractional multi-order heat conduction model of the human head

The subject of the article is devoted to the development of a matrix collocation technique based upon the combination of the fractional-order shifted Vieta–Lucas functions (FSVLFs) and the quasilinearization method (QLM) for the numerical evaluation of the fractional multi-order heat conduction model related to the human head with singularity and nonlinearity. The fractional operators are adopted in accordance with the Liouville–Caputo derivative. The quasilinearization method (QLM) is first utilized in order to defeat the inherent nonlinearity of the problem, which is converted to a family of linearized subequations. Afterward, we use the FSVLFs along with a set of collocation nodes as the zeros of these functions to reach a linear algebraic system of equations at each iteration. In the weighted \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2$$\end{document}L2 norm, the convergence analysis of the FSVLFs series solution is established. We especially assert that the expansion series form of FSVLFs is convergent in the infinity norm with order \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {O}(\frac{1}{K^3})}$$\end{document}O(1K3), where K represents the number of FSVLFs used in approximating the unknown solution. Diverse computational experiments by running the presented combined QLM-FSVLFs are conducted using various fractional orders and nonlinearity parameters. The outcomes indicate that the QLM-FSVLFs produces efficient approximate solutions to the underlying model with high-order accuracy, especially near the singular point. Furthermore, the methodology of residual error functions is employed to measure the accuracy of the proposed hybrid algorithm. Comparisons with existing numerical models show the superiority of QLM-FSVLFs, which also is straightforward in implementation.

The study of the heat conduction associated to the human head subject to various environmental temperatures has been attracted by several research scholars over the past decades.From the theoretical point of views, the temperature distribution of the head based on physically realistic assumptions were investigated in 1,2 .This process was modeled as a two-points boundary value problem (BVP) with nonlinearity and singularity, see also 3 .In a simplified form, the governing equation in 1-D steady state and in the spherically symmetrical form is described as follows subjected to the boundary conditions given as follows www.nature.com/scientificreports/Here, T(r) stands for the absolute temperature depending on r as the radial distance measured from the center to the periphery of the head.Moreover, by K we denote the (average) thermal conductivity, β represents the heat exchange coefficient, δ shows the thermogenesis heat production factor, and α is the metabolic thermogenesis slope parameter.Finally, T a is denoted the ambient temperature.Next, we introduce the following dimensionless variables given by By utilizing the former relations in (1.1), we arrive at the integer-order model Here, the new parameter Bi denotes the Biot number, the parameter represents the thermogenesis heat production factor, and the metabolic thermogenesis slope parameter denoted by m .For more information about the preceding equations (1.1) and (1.3) we refer to 1-3 .Generally, one can not find an analytical or a closed-form solution to (1.3) due to the (exponential) nonlinearity term and two different parameters m, in (1.3).As a consequence, both semi-analytical and numerical solving approaches become a preferred alternative.Various strategies have been developed in the literature for the BVP (1.3) or the related type equations.In 3 , the varational technique proposed to find solutions of (1.3) accurately.By utilizing the maximum principles, an analytical pointwise lower and upper bounds for the exact solutions derived in 4 .The author of 5 provided a semi-analytical and non-perturbative solution to (1.3).The stochastic scheme utilizing artificial neural networks investigated in 6 .Some other computational techniques have been developed for (1.3) such as the varational iteration scheme 7 , the sinc-Galerkin procedure 8 , the reproducing kernel method 9 , the modified homotopy analysis approach 10 , and the compact finite difference scheme 11 .
In contrast to the traditional definition of the derivative, various concepts of fractional-order derivatives have been defined in the literature.Indeed, a fractional-order system can provide memory effects and produces solutions, which are much closer to the nature of real-life phenomena, see 12,13 .For more information and applications of various fractional models we refer to [14][15][16][17][18][19][20] .In this work, we are interested to interpret the derivatives involved in (1.3) in the sense of Liouville-Caputo derivative.To be more precise, the following class of fractional multi-order BVP with singularity and nonlinearity is considered Here, the derivative operators LC D γ r , LC D σ r are taken in the sense of Liouville-Caputo of orders γ ∈ (1, 2] and σ ∈ (0, 1] respectively.Moreover, the parameter A is a real constant.Also, two model parameters m, are previ- ously defined in (1.3).It should be emphasized that in addition to the existing nonlinearlty and singularity, the fractional orders γ , σ make the proposed model more complicated than (1.3).These hinder us to find a closed- form solution to this class of equations.Therefore, proposing numerical solutions are preferred to deal with (1.4).To the best in our knowledge, the only method developed is the polynomial least-squares scheme 21 with A = 2 .On the other hand, the human head model using the He's fractal derivative was considered in 22 .The resulting fractional model solved by the Taylor expansion series together with two-scale transform technique.
The aim of this research study is to design not only an efficient but also an accurate procedure based on a novel family of polynomial functions to solve (1.4) numerically.Our strategy is relied on a matrix collocation technique utilizing a generalized version or a fractional-order counterpart of the basis functions.The successful applications of the collocation methods with exponential convergence properties and using various bases have been proven in solving many important and practical problems.Among others, we refer to collocation matrix procedures based on Bessel, Chebyshev, Jacobi, and Legendre functions [23][24][25][26][27][28] , to name a few.The bases will be used here are the Vieta-Lucas polynomials 29 .The shifted form of this set of orthogonal functions has been recently defined in the literature, see cf. 30,31 .However, in the present work we are pioneering in defining and utilizing the fractional-order version of shifted Vieta-Lucas functions (SVLFs) to acquire the approximate solution of (1.4) accurately.On the other hand, we employ the quasilinearization method (QLM) to the model equation (1.4) to have efficiency in our proposed algorithm.The second advantage of the QLM is that we have no longer any nonlinearity term in the model.In other words, we have to solve a set of quasilinear equations in an iterative fashion.Besides the aforementioned benefits of the presented QLM-FSVLFs approach, we also establish the uniformly convergence of the FSVLFs in the L ∞ norm.The same convergence property is proved with regard to the L 2 norm.
The outline of the current study is considered next.Some concepts of fractional calculus is provided in "Fractional derivative: the Liouville-Caputo operator" section.The definitions of the Vieta-Lucas functions and their shifted and generalized forms are done in "The shifted Vieta-Lucas functions: their generalized form and error estimate" section.The proofs of the convergence of FSVLFs in both weighted L 2 and L ∞ are then given.In "The combined QLM-FSVLFs procedure" section, we give the details of the main proposed combined QLM-FSVLFs algorithm.The test of accuracy of the presented collocation procedure is accomplished by using the definition of the residual error functions at the end of this section."Numerical simulations" section is devoted to the computational results in order to verify the efficacy as well as the accuracy of the presented algorithm.Comparison are also made with numerical calculations of available computational algorithms.The last "Conclusions" section contains a conclusion of the study. (1.3) r ∈ (0, 1).

Fractional derivative: the Liouville-Caputo operator
Let us review some fundamental concepts from fractional calculus theory.They are essential in our subsequent discussions in the following sections.For a detailed description we refer the interested readers to the monograph and standard textbook 12,13 .First, note that by Ŵ(•) we denote the well-known Gamma function.We say that a real-valued function p(r), r > 0 belongs to the space C η , η ∈ R if there exits a function k(r) ∈ C([0, ∞)) and a real number τ > η such that we have p(r) = r τ k(r) .Furthermore, we call p(r) ∈ C m η if and only if p (m) (r) ∈ C η for a m ∈ N. Now, we are in a position to state the definition of the integral operator in the sense of Riemann-Liouville (RL).Let suppose that p(r) ∈ C η , η > −1 .The RL integral of function p(r) of order γ > 0 is defined by Next, we give the definition of Liouville-Caputo (LC) derivative.One further remarks that the fractional LC operator LC D γ r has the linear property.Besides, let C be a real constant number, then one has The following property computes the fractional LC derivative of p(r) = r w .Thus, we have Here, N 0 := N ∪ {0} .Moreover, ⌈ξ ⌉ and ⌊ξ ⌋ stand for the ceiling and floor functions.While the former outputs the smallest integer number greater or equal than ξ , the latter function represents the largest integer number less or equal than ξ.

The shifted Vieta-Lucas functions: their generalized form and error estimate
We begin this section by introducing the definition of shifted Vieta-Lucas functions (SVLFs).These functions were originally introduced in 29 and their shifted form recently utilized in some papers 30,31 .In the next stage, we define the generalized version of SVLFs.In the last part, we examine the error estimate as well as the convergence properties of SVLFs in a rigorous way.

The SVLFs and their generalization
The Vieta-Lucas functions on the interval [−2, 2] are defined via the following relation Obviously, we have V 0 (s) = 2 .Also, a little manipulation shows that V 1 (s) = s .By making use of s = 4r − 2 , we get the shifted form of these functions for r ∈ [0, 1] .Let us denote the new version by V ⋄ k (r) , which is equal to V k (s) .By setting V ⋄ 0 (r) = 2 and V ⋄ 1 (r) = −2 + 4r , the following recursive formula gives us the remaining SVLFs as From (3.2) we get the next two terms of these functions as Thus, at r = 0, 1 as the special values we arrive at The proof of the former relations can be easily done by induction on k on the recursive formula (3.2).
Note that the SVLFs solves the following differential equations.In the Sturm-Liouville representation, we can write it as where w(r) := r − r 2 − 1 2 is known as the weight function.From relation (3.3), one can prove that these functions are orthogonal with regard to w(r) on (0, 1) .The orthogonality condition is The explicit analytical form of the SVLFs is given by By setting V ⋄ k (r) = 0 , we can determine the zeros or roots of SVLFs.The zeros of V ⋄ k (r) of degree k ∈ N are all real and distinct on the open interval (0, 1) given by where s k p := 2 cos are the zeros of normal Vieta-Lucas functions V k (s) defined in (3.1).Following the previous works 23,32 , the main goal is to employ generalized version of SVLFs to get more accurate results in numerical computations.The next definition will do this task.Definition 3.1 The fractional-order SVLFs (FSVLFs) are denoted by V α k (r) on [0, 1] of degree k.They are defined through making use of the change of variable r → r α , α > 0 as With the aid of the aforementioned change of variable in (3.5), the following analytical form is obtained By a little calculation one can show that the family {V α k (r)} ∞ k=0 constitutes a orthogonal set and the related weight function is w α (r) = r α−1 √ r α −r 2α for r ∈ (0, 1) .Thus, we get The roots of FSVLFs can be utilized as a candidate for the collocation points in our matrix collocation algorithm.Using the fact that the SVLFs of degree r can be expressed as its zeros shown in (3.6), we have Therefore, the roots of V α k (r) are defined by r 1/α p for p = 1, 2, . . ., k .Thus, we have proved the following fact: Lemma 3.2 The roots of the FSVLFs, V α k (r) of degree k, are within (0, 1) and given by for p = 1, 2, . . ., k.

Analysis of error and proof of convergence
Next, we are interested to investigate the behavior FSVLFs in an rigorous way.To begin, we recall that a square integrable function g(r) may be written as a linear combination of FSVLFs.This implies that To acquire the unknown coefficients β k , we can employ the orthogonality relation (3.9) to reach where η 0 = 4 and η k = 2 for k > 0 .Let us further denote by h 2,w the weighted L 2,w norm on [0, 1] with regard to weight function w α (r) .To establish that the series solution (3.11) is convergent in the L ∞ norm, one requires (3.4) (3.9) (3.10) (3.12) www.nature.com/scientificreports/ to prove that the coefficient β k are decreasing functions of k .To this end, we first derive an upper bound for β k in (3.11).Let us emphasize that the following results are valid for α = 1 .For a proof related to α = 1 we refer to 33 .
Then, an upper bound for β k in (3.12) is given as Proof By making use of change of variable r = 1 2 (1 + cos s) or r = cos 2 s 2 =: p(s) in (3.12) to render Integration by parts (twice) yields where By using the upper bound for the second derivative we arrive at It is now remained to find the value of the integral term in (3.16).By changing of variables v = (k ± 1)s , we get for an odd k > 1 as . Consequently, we arrive at By inserting (3.17) into (3.16) and using the fact that η k = 2 for k ≥ 1 , the proof is finished.
In practice we have to utilize a cutted series representation to approximate g(r) instead of the infinite series form (3.11).By truncating (3.11) up to its first (K + 1) terms, we have Let assume that two consecutive approximations to g(r) will be denoted by g K (r) and g K+1 (r) .Our goal is to measure the difference between them as In the next result, we derive an upper bound for the error e K in the weighted L 2 norm.

Theorem 3.4 Let suppose the assumptions of Theorem 3.3 are valid. Then, we get an error estimate as
Proof We first exploit the definition of error (3.19) followed by using (3.18) to get (3.15) (3.17) www.nature.com/scientificreports/Applying the orthogonality condition (3.9) reveals that �V α K+1 (r)� 2,w = √ π/(2α) for K ≥ 1 .Hence, the result of Theorem 3.3 with α = 1 gives us The error between the cutted (finite) series form g K (r) in (3.18) and the infinite expansion series (3.11) will be defined next by We refer to this error as the global error.We are interested to estimate an upper bound for E K (r) in both L 2,w and L ∞ norms.The first upper bound is obtained in the weighted L 2,w ([0, 1]) norm.The inequality (3.13) given in Theorem 3.5 is now applied to the foregoing formula.This implies by using α = 1 that If one uses the Integral Test 34 , the immediate conclusion is We then insert the last inequality into (3.21).Ultimately, we take the square root to get the desire assertion.

Remark 3.6
It should be emphasized that our error bound is of order O(K −7/2 ) .This upper bound is an improvement over the existing upper bounds in the literature 31,35 .In fact, they obtained the error bound as In particular, we will prove that the norm of the error �E K � ∞ = O(K −3 ) .This shows that the expansion series solution g K (r) converges uniformly as K → ∞.
Still we are interested to estimate the global error (3.20) in the L ∞ norm.Theorem 3.7 Suppose that the assumptions of Theorem 3.3 satisfied.Then, the (global) error

has the following bound
Proof By employing the triangle inequality we get . Vol

The combined QLM-FSVLFs procedure
We note that the technique of matrix collocation with the aid of FSVLFs can be utilized to solve the nonlinear human head model (1.4) directly.However, the main deficiency of applying the direct method is that the resultant system of equation has nonlinearity.This implies that the spent time to solve the nonlinear system increases as a function of K, the number of FSVLFs bases, in order to attain an acceptable degree of accuracy.To overcome the inherit nonlinearity, we propose to apply the technique of quasilinearization to the original model (1.4).It transforms the model under consideration into a sequence of linearized subequations.In accordance to theory of the Newton method by using a rough first approximation, we get a quadratic convergence to the solution of the heat conduction model related human head in model problem (1.4).
Let us describe the basic facts related to the quasilinearization method (QLM).After converting the nonlinear model of human head (1.4) into a set of quasilinear submodels, the direct FSVLFs matrix collocation strategy is applied to every linearized models.In what follows, we refer to this hybrid technique as the QLM-FSVLFs approach.For detailed descriptions and also the recently published papers on the applications of QLM, see cf. [37][38][39][40] .
To begin with, we rewrite the nonlinear human head model (1.4) as Here, L γ , denotes the linear operator and G is a nonlinear function.They defined by To proceed, we suppose that H 0 (r) denotes the rough first approximation to H(r) as the true solution of (4.1).Now, the process of QLM for (4.1) can be given as Along with the preceding equations, the same initial conditions as in (1.4) are given by Doing some straightforward computations, the application of QLM yields the next expression for (4.1) For the sake of abbreviation, we make use of the following notations for ℓ = 0, 1, . . . to write the last equation (4.3) in a compact expression as The solution of the quasilinear multi-order BVPs (4.4) can now be determined with the help of the FSVLFs matrix collocation strategy.Using the same expansion series (3.18), the approximate solution of (4.4) is taken as a cutted series with (K + 1) bases as  K,α (r) .To do so, we have to calculate the corresponding high-order derivatives of the vector R R R α K (r) .Owing to the boundary conditions (4.5), we need to compute only the first-order derivative of the unknown solution.Additionally, the other requirements are the computations of the fractional derivatives of order γ and σ as they are appeared in (1.4).In this respect, we utilize two properties (2.1) and (2.2) as basic tools.According to the previous works 32 (Algorithm 3.1) or 23 (Algorithm 4.1), we are able to calculate the γ -and σ-derivatives of R R R α K (r) efficiently and straightforwardly.The complexity of the proposed algorithm is bounded by O(K + 1) .Henceforth, the vector representation forms of the derivatives are denoted by Once we have the derivatives of R R R α K (r) at hand, we can easily deduce the followings: (a) The first-order derivative of h K,α (r) in (4.8) is given by To continue, we utilize the zeros of FSVLFs as defined in (3.10).Note that we have (K + 1) unknown coef- ficients in (4.5).Thus, the zeros of V α K+1 (r) will be utilized as these are labeled as r 0 , r 1 , . . ., r K .We are now in a position to place these nodes into the vector formats of the unknown solutions and its fractional-order derivatives.So, we introduce the following matrix forms We thus reach at the next assertion: Lemma 4. 5 In the vector formats, the approximate solution h (ℓ+1) K,α (r) and its fractional ψ-derivatives LC D ψ r h (ℓ+1) K,α (r) for ψ = σ , γ computed at the zeros of FSVLFs can be represented as where the vectors of R R R and R R R ψ have already defined by (4.11) and Let us come back to the quasilinear form (4.4).We then insert the roots of FSVLFs into (4.4).By employing relations (4.13), one arrives at the next matrix form where we have used the following notations We next employ two relations (4.12).At each iteration ℓ , we get the fundamental matrix equation given by In a more convenient form we have where for ℓ = 0, 1, . . . .One still requires to incorporate the given boundary conditions (4.2) into the matrix format.Hence, they will be added to the fundamental matrix relation (4.16).Firstly, we consider the boundary condition h ′ ℓ+1 (0) = 0 .It is sufficient to employ the relation (4.9).If we let r approaching to zero, we have the matrix format For the endpoint boundary condition h ′ ℓ+1 (1) + Bi h ℓ+1 (1) = Bi , we must use both relations (4.8) and (4.9).By approaching r → 1 we obtain the matrix format (4.9) d dr h www.nature.com/scientificreports/By using these two row matrices [ X X X 0,ℓ ; 0] and [ X X X 1,ℓ ; Bi] we will substitute two rows of the fundamental matrix equation [X X X ℓ ; � � � ℓ ] in (4.16).We denote the resulting modified system as After solving (4.17), we get the unknown coefficients β (ℓ) k for k = 0, 1, . . ., K at each iteration ℓ = 1, 2, . . . in (4.5).We remark that usually taking r = 5 is enough to reach the desired accuracy when we solve the system (4.17) through utilizing the QLM-FSVLFs technique.

The methodology of residual error functions (REFs)
Practically, the analytical or exact solutions to (1.4) are not tractable.Particularly, if we dealing with the fractional orders 1 < γ ≤ 2 and 0 < σ ≤ 1 .Therefore, the approach of REFs is employed to calculate the accuracy of the developed QLM-FSVLFs.For this purpose, we insert the obtained approximate solution (4.5) into (1.4).Thus, the REFs given by for ℓ = 0, 1, 2, . . . .Next, the maximum values the REFs (for a fixed ℓ ) is calculated by defining Finally, we calculate the notion of numerical order of convergence ( ord ∞ K ).It is defined by the ratio

Numerical simulations
The goal is here to describe the performance of the developed matrix collocation approach based on the FSVLFs to the human head model (1.4) with singularity and multi-order fractional For this purpose, several numerical computations based on MATLAB software version 2021a have been provided to show the basic fundamental results obtained in this manuscript.Diverse values of parameters , m as well as different fractional orders γ , σ are considered to illustrate the utility of the QLM-FSVLFs technique.Below, the parameter ℓ in the QLM is set as ℓ = 5 in the practical computations.We further use the initial approximation H 0 (r) as the zero function.
The integer-orders γ = 2 and σ = 1 As a starting point, we first consider the integer-order derivatives.Note that we always use α = 1 when the deriva- tives are both integers.To continue, we take the next parameters for the human head model (1.4) as By utilizing the QLM-FSVLFs procedure with K = 6 we obtain The above approximate solution is depicted in Fig. 1.To validate our finding, we compare our outcomes with those reported by the polynomial least-squares method (PLSM) 21 with the same number of basis functions; K = 6 .The reported approximation is as follows A simple comparison indicates that the results of both methods are very close together.This fact can be further justified by computing the difference between two approximate solutions via |h 6,1 (r) − ũ(r)| .This error is shown in Fig. 1, on the right panel.We also visualize the REFs R We next utilize larger values of K = 8, 10, 12, and K = 14 in the computational experiments to justify that the presented QLM-FSVLFs approach converges numerically toward to the actual solution.The achieved REFs related to these values are presented in Fig. 2. It can be readily seen that the REFs is a decreasing function of K as increased.
We further validate the computational outcomes gotten by the above aforementioned parameters.To do so, we compare the outcomes with those reported by the well-established existing numerical models as observed in Tables 1 and 2. The first approach is the PLSM 21 with polynomial of degree six while the second one is the stochastic based on an unsupervised artificial neural networks optimized with genetic algorithms (GAs), active-set technique (AST), and their hybrids named GA-AST from 6 with ten neurons.Our presented results are obtained by utilizing K = 6, 8 , and K = 10 .For justification, the related achieved REFs are also tabulated in these tables.One can clearly infer that our results with K = 6 are comparable with the outcomes of the PLSM in terms of 6,1 (r) = −0.00002031222803r 6 + 0.00000919680423 r 5 − 0.0008265202274 r 4 + 0.000004263676977 r 3 − 0.05220576355 r 2 + 1.160819842.www.nature.com/scientificreports/accuracy.However, our results are more accurate using K = 8 with those errors reported by the ANN-based procedures in Table 2.
The behavior of the obtained errors E ∞ in the QLM-FSVLFs are also examined using the above model parameters as seen in Table 3.The corresponding numerical order of convergence, i.e., ord ∞ R , are reported alongside in Table 3.Here, we set K = 2, 4, 8, 16 and the relation (4.20) is utilized.In the case of m = 1 and for each K , the (spent) CPU times measured in seconds are further presented in Table 3.One should note that the needed CPU time is just the time required to solve the modified fundamental system of equations [ X X X ℓ ; � � � ℓ ] in (4.17).The results presented in Table 3 indicate that the developed QLM-FSVLFs with a linear cost produces high-order accurate outcomes.
Beside the results shown in Table 3, we next investigate the impact of employing different model parameters , A , and Bi while the others are assumed to be fixed.To do so, we first fix m = 1 , A = 2 , and Bi = 1 .Figure 3 presents the numerical solutions obtained by various = 0.25, 0.5, 1, 2 , and = 5 .The number of bases used here is K = 10 .The achieved REFs associated to these approximate solutions are further visualized in Fig. 3.It can be clearly seen that the REFs are increasing as is increasing.Next, we fix m, = 1 , A = 2 and vary Bi = 0.25, 0.5, 1, 2, 5 .The results are shown in Fig. 4. Here, it can be also seen that the REFs are increasing if we increasing the value of the Biot number Bi from 0.25 to 5. In the last experiment, we visualize the numerical solutions for different values of the coefficient parameter A = 0.25, 0.5, 1, 2, 5 .In these cases, the other parameters are assumed to be fixed as m, , Bi = 1 .Figure 5 shows the aforementioned visualizations.Contrary to the other parameters, the REFs will be decreased as A increases.
Lastly, for the worst case scenarios, i.e. = 5 , Bi = 5 , and A = 0.25 we compute the maximum absolute errors E ∞ as given by (4.19) and the associated numerical order of convergence ord ∞ K defined by (4.20).These outcomes are reported in Table 4.One can easily conclude the high-order accuracy of the developed QLM-FSVLFs technique.

The fractional-order derivatives
Here and in the following part, let us consider the fractional-order derivatives γ and σ in the human head model (1.4).Furthermore, the effect of simultaneous use of different parameters m, , Bi, A , and γ , σ on the calculated approximate solutions will be investigated.In particular, the influence of the local parameter α is of important for us.We will present these impacts through figures and tables.For the fractional multi-order model problem (1.4), let us first consider to the next parameters as they recently utilized by 21         Firstly, let us consider K = 6 and α = 1 .By utilizing the QLM-FSVLFs the approximate solution using these parameters is given by Similar to the integer-order cases, let us mention the numerical results reported by the polynomial least-squares method (PLSM) 21 with the same number of polynomial basis functions.This approximation was given as These two approximate solutions are plotted in Fig. 6, on the left picture.To be more precise, we also show the graphical representation of the error between our numerical solution and that solution obtained by PLSM.This error is depicted in Fig. 6, at the right plot.
Let us compare the errors obtained by using the three different choices of the (local) parameter α equal to 1, σ , γ .Here, we have γ = 1.7 and σ = 0.7 .Motivated by the previous works 23,32 , we also here show that the achieved residual error becomes smallest in magnitude as one selects α = γ .By taking K = 6 , these results are presented in Fig. 7. To see choosing the α = γ leads to the best possible result, we further plot the E ∞ obtained by the QLM-FSVLFs technique versus the parameter α as depicted in Fig. 8.These results are calculated with K = 10 and α ∈ [0.05, 3] with step size 0.05.Indeed, the approximate solution using α = γ = 1.7 and σ = 0.7 with K = 6 is as follows A comparison is also done in Table 5 when we have the fractional-orders γ = 1.7 and σ = 0.7 .In these experiments, we utilize K = 6 , m, , Bi = 1 , and A = 2 as before.For validation, the outcomes of the previously well-established scheme, namely the PLSM 21 are shown in this table.One can see that by employing the parameter α to be as the fractional order γ leads to more accurate results.Our results using this value of α are considerably more accurate than those reported by PLSM.
.8809 × 10 −3 5.5447 7.3924 × 10 −3 www.nature.com/scientificreports/ that our presented spectral matrix collocation technique with less computational efforts produces more accurate outcomes.Ultimately, we fix the fractional-order parameters given by Then, we consider the impacts of diverse parameters m, , Bi on the approximate solutions of the model equa- tion (1.4) through the QLM-FSVLFs method.To this end, we set K = 10 and use α = γ , which is 1.75.Also, we take A = 2 in all subsequent experiments.First, we vary = 0.25, 0.5, 1, 2, 5 and fix other parameters as m, Bi = 1 .The results of approximations are plotted on the left part of Fig. 10.Furthermore, on the right part, the results of REEFs are presented.As in the integer-order cases, by increasing the nonlinearity parameter , the achieved errors get increased.In the next experiments, we will consider the effect of the nonlinearity parameter m.In this respect, we vary m = 0.25, 0.5, 1, 2, 5 .The outcomes of numerical solutions along with the associated REFs are shown in Fig. 11.Obviously, the magnitude of the achieved errors will be decreased if the factor of nonlinearity m is increased in the model of human head.
In the final stage, let us investigate the effects of using the boundary condition parameter Bi on the calculations of QLM-GFSVLFs.The L ∞ error norms (4.19) together with the corresponding numerical order of convergences calculated via (4.20) are further shown in Table 8 for the values of m = 0.25, 0.5, 1, 2, 5 .It can be evidently seen that the developed QLM-FSVLFs technique produces high-order accurate outcomes when applied to the model (1.4).

Conclusions
A numerical spectral collocation procedure relied on the (novel) fractional-order shifted Vieta-Lucas functions (FSVLFs) designed to get the polynomial solutions to a class of nonlinear fractional-order two-point boundary value problems with singularity arising in modeling of heat conduction of the human head.Both fractional derivative operators are interpreted in the sense of Liouville-Caputo's derivative.To conquer the inherit nonlinearity of the model problem, the quasilinearization method (QLM) is utilized to reach at a family of linearized γ = 1.75, σ = 0.75.subequations.Hence, we solved this family of equations through using the spectral matrix collocation scheme relied on the FSVLFs.Theoretically, we proved that the related series solution of FSVLFs is convergent in the weighted L 2 norm.In particular, we showed that the series of FSVLFs is uniformly convergent of order O(K −3 ) , where K is the number of FSVLFs used in the approximation.Simulation results are provided to verify our theoretical findings and illustrate the effectiveness of the presented QLM-FSVLFs scheme.Numerical results are well tabulated through tables against the various model parameters m, , Bi, A as well as fractional orders γ and σ .We have compared our outcomes and plots with those reported by the PLSM 21 and there is a good level of agreement in our results.The reported rate of convergence shows that the QLM-FSVLFs method converges exponentially as the number of bases increases.The proposed QLM-FSVLFs technique with high-order accuracy can be easily generalized to tackle various multi-order multi-dimensional fractional model problems with applications in various disciplines of science and engineering.

Definition 2 . 1
Let m − 1 < γ < m , m ∈ N and assume that p(r) ∈ C m −1 .The LC fractional derivative of p(r) is defined by where D = d dr .

Theorem 3 . 5
Under the assumptions of Theorem 3.3, the (global) error term E K (r) in the L 2,w ([0, 1]) norm is bounded as Proof In accordance to the definitions (3.11) and (3.18) we get The orthogonality condition (3.9) yields

Lemma 4 . 1 7 ) 4 . 2
(r) in(4.5) in the formThe next Lemma provides a decomposition for V V V α K .Its proof is easily obtainable by considering the relation (3.8) in Definition (3.1).The vector of FSVLFs is represented by where the vector of monomials is defined by and N N N K as an upper triangular matrix has the following structure The next corollary is obtained by combining two former relations (4.6) and (4.Corollary The approximate solution h(ℓ+1) K,α (r) in (4.5) is rewritten as Remark 4.3 It should be stressed that the matrix representation (4.8) is beneficial when we need to have the high-order derivatives of h (ℓ+1) https://doi.org/10.1038/s41598-024-53822-6www.nature.com/scientificreports/(b) The fractional-order derivatives of LC D (r) in (4.8) for ψ = γ , σ are computed as (r) achieved via (4.18) as depicted in the former figure.
34 now use the fact that the Vieta-Lucas functions can be expressible in terms of Cheyshev function of the first kind29Also, we have that36|T r (s)| ≤ 1 for all s ∈ [−1, 1] and for all r.Therefore, by changing the variable 4r α − 2 = 2 cos s we find that |V α k (r)| ≤ 2 for all r ∈ [0, 1] .Using the upper bound derived in (3.13) we get , we apply the Integral Test34.It follows that Now, the last identity is placed into the inequality (3.23).If we take the supremum over all r ∈ [0, 1] , we get the desired assertion.