Computational power is a challenge when it comes to the high-fidelity modeling of nuclear reactors. Detailed simulations of reactor physics involve complex calculations that require significant computing resources, which can be time-consuming and expensive. Reduced Order Modeling (ROM) allows for an approximation of a complex model by only capturing the essential features, thereby reducing the computational load. A reduced order model provides computationally efficient approximations of a system, but it requires still many evaluations of a high-fidelity model to capture all the dynamics. Using the adaptive sparse grid can reduce the number of evaluations needed, though the construction of the reduced order model is still computationally intensive. The aim is to minimize the computational workload involved in constructing a reduced-order model during the offline phase. This is achieved by decreasing the number of high-fidelity model evaluations necessary for building the reduced order model while maintaining accurate results. To this end, the existing adaptive proper orthogonal decomposition algorithm is enhanced by employing multi-fidelity techniques. Multi-fidelity methods aim to combine large amount of low-fidelity data with a limited amount of high-fidelity data to compute accurate, yet computationally inexpensive approximations. Two novel multi-fidelity reduced order model methods based on proper orthogonal decomposition are proposed; Filtered Bi-Fidelity Adaptive Proper Orthogonal Decomposition (FB-POD) algorithm and Adapted Bi-Fidelity Proper Orthogonal Decomposition (AB-POD). These models are evaluated on two different test cases, and the balance between the accuracy of each multi-fidelity ROM and the computational cost, measured by the number of high-fidelity evaluations, is investigated. In specific cases, the proposed methods significantly reduce the number of high-fidelity evaluations compared to the single high-fidelity ROM, while yielding comparable accuracy.

The growing amount of plastic pollution in the oceans is worldwide a huge problem. There is an increasing amount of projects that try to clean up this pollution. To remove the pollution efficiently, it is useful to know how plastic litter moves, degrades and sinks. Models are being developed that can describe what happens with plastic in the ocean. To model the behaviour of plastic in the ocean a lot of parameters need to be taken into account. The transport of plastic particles is largely determined by processes like the wind and the current but their movement also has a random character due to dispersion. Due to this random movement, the future position of a particle can be described by a probability distribution. To find that probability distribution, the Kolmogorov forward equation (or Fokker-Planck equation), in this context the same as the advection-diffusion equation, needs to be solved. It is also possible that the plastic particles sink and settle on the bottom of the ocean. One goal of this project was to incorporate the settling of particles as an extra term into the advection-diffusion equation. Forward models are used to predict where particles will go to. It can also be useful to have reverse time models that can describe where particles had their origin. To make reverse time models, it is necessary to solve a reverse time advection-diffusion equation. The main goal of this project was to derive and solve the reverse time advection-diffusion equation that includes the settling of particles. This is done by finding the adjoint of the Kolmogorov forward equation and the settling term.

This thesis investigates the tissue sparing effect of FLASH (>40 Gy/s) radiation, as opposed to CONV (conventional, dose rates typically between 0.01­0.1 Gy/s) radiation. We irradiated zebrafish embryos (4 days past fertilisation) with 116 MeV protons. The aim was (1) to measure the effect, and (2) if the effect were significant, see whether it depended on the oxygen concentration in the tissue, as the oxygen depletion hypothesis (a popular theory on the underlying mechanics of the FLASH effect) predicts. We irradiated embryos with either FLASH or CONV, where a possible FLASH effect would reduce toxicity of the FLASH radiation. We did the same for zebrafish which were deliberately put in a hypoxic condition prior to irradiation. In that case, the depletion hypothesis would predict that the difference between FLASH and CONV disappears. Our biomarkers for radiobiological damage were the survival rate and γ­H2AX foci formation. In our experimental conditions, the radiation effect on the survival rate was eclipsed by other factors which could not be isolated. We did confirm the possibility of using γ­H2AX foci formation as a marker for radiobiological damage in full­body irradiated zebrafish embryos. There were individual samples that showed clear and localised specific γ­H2AX signal, but these were too scarce and the signal was too inconsistent across samples to gather meaningful statistics. This was most often caused by limited antibody penetration in the embryo. We were therefore unable to draw conclusions about the FLASH effect. Better and more consistent antibody penetration, e.g. by longer digestion in collagenase before antibody staining, could change this in the future. We further custom­built and validated a hypoxic aquarium to produce hypoxic zebrafish tissue, as well as a computational model of the irradiation setup to simulate the dose distribution in the zebrafish container. We found the dose distribution to be sufficiently homogeneous for our experiment, at least 91.47% uniform for CONV and 90.72% uniform for FLASH.

The delineation of the regions of interest (ROIs) plays an important role in the radiotherapy workflow. The ROIs include the gross tumor volume (GTV), the clinical target volume (CTV), the planning target volume (PTV) and organs at risk (OARs). There are however uncertainties related to the delineation of the ROIs due to inter­observer variability between radiation oncologists caused by e.g., a lack of consensus on anatomic definition or a lack of contrast in medical images. The largest inter­observer variability has been found in esophageal cancers, head and neck cancers, and lung cancers. To pre­ vent inter­observer variability, auto­contouring software can be used to delineate the ROIs. There are however also uncertainties in the delineations made by auto­contouring software. The purpose of this study was to perform an accurate evaluation of the dosimetric effect of delineation uncertainties. To do so, first a method to characterize the delineation uncertainties had to be found. The delineations uncertainties were characterized with principal component analysis (PCA). By performing PCA on a set of delineations, the eigenmodes which represent the variations between the delineations can be found and new random delineations based on the eigenmodes can be formed. To determine the dosimetric impact of the delineation uncertainties, two cases were investigated: (1) the dosimetric impact of delineation uncertainties for a fixed dose distribution; and (2) the dosimetric impact of delin­ eation uncertainties for a dose distribution reoptimized for every possible realization of a delineation. For the fixed dose distribution, Polynomial Chaos Expansion (PCE) was used as a meta­model for the dose volume histogram (DVH) of the target, and for the reoptimized dose distribution PCE served as a meta­model for the total dose distribution and the DVHs of the target and other ROIs. The delineation uncertainties of two data sets were analyzed: (1) 12 manual delineations of the GTV of a hepatocellular carcinoma patient; and (2) 90 auto­contours of the CTV and brainstem of a head and neck patient. The variation of the manual delineations of the GTV could be accurately described by 5 eigenmodes, while 45 eigenmodes were needed to describe the variations of the brainstem auto­ contours. The delineation uncertainties in the auto­contours of the CTV could not be characterized by PCA due to the shape of the CTV. The dosimetric effect of the delineation uncertainties of the manual delineations for a fixed dose distri­ bution was investigated for both an intensity modulated proton therapy (IMPT) plan and a volumetric­ modulated arc therapy (VMAT) plan for 10,000 random delineations of the CTV. For this patient the CTV was equal to the GTV. In the VMAT plan, the CTV received sufficient dose for all delineations, but the PTV was underdosed in 17.1% of the delineations. In the IMPT plan, the CTV was underdosed in 69.0% of the delineations. This percentage of underdosed CTV delineations in the IMPT plan was reduced when the plan was made robustly. The results for the fixed IMPT dose distribution should however be verified with a more accurate PCE model. A proof of principle for the PCE as a meta­model for the reoptimized dose distribution and DVHs of the target and other ROIs was shown. However, a more accurate PCE model with a higher polynomial order would be needed to analyze the effects of the reoptimization on the dose delivered to the ROIs. It has been shown that the dosimetric impact of delineation uncertainties can be modelled using PCE. This is a first step towards systematically and quantitatively taking into account delineation uncertainties in radiotherapy treatment planning. In future research the analysis of the dosimetric impact of the un­ certainties can e.g., be used in adaptive radiotherapy in which auto­contouring is used. By knowing the dosimetric impact of the uncertainty, treatment plans could be optimized such that a structure receives its target dose with a certain probability or critical areas where the dosimetric impact of the delineation uncertainties is large could be flagged such that these areas are checked before treatment.

Pencil beam scanning is becoming a more common treatment modality. However, its ability to deal with moving targets is known to be limited, as beam motion and target motion can reinforce each other, deteriorating the planned dose distribution in what is called the interplay effect. Literature concerning breathing motion usually investigates regular patterns. This work aims to investigate the magnitude of the interplay effect when considering irregular breathing signals. In silico calculations of dose distributions were made in the treatment planning system RayStation (version 7.99), using an XCAT phantom with 50 CT phases to model moving patient anatomy. An interplay calculator was included in RayStation, allowing calculation of disturbed doses based on a treatment plan and an irradiation time model for a proton therapy accelerator. The target investigated was a spherical liver tumour with 5cm diameter, irradiated with two beams delivering a prescription dose of 63 Gy. Plans without and with 5x layered repainting were created. Clinically realistic regular breathing patterns were generated to establish a baseline, after which irregularities were introduced. The basic form for all patterns was a sin^4 signal, with regular signal amplitudes ranging from 6 to 18 mm, period ranging from 3 to 4 s and phase between 0 and 2π rad. Considered irregularities were baseline shifts up to 34 mm, changing amplitudes between 6 and 18 mm, changing periods between 1.6 and 5.2 s and combinations. Evaluation was done by looking at dose homogeneity HI_5 and the fraction of the CTV volume that received a dose outside of the clinical limits of 95% and 107%, V_107/95. For the regular patterns, both a systematic and a randomised analysis were carried out. For irregular patterns, only a systematic analysis was carried out. The mean HI_5 was found to be 31% for regular patterns; the means of all irregular patterns stay below this, even though the size of the irregularities for some breathing patterns was very large. The mean V_107/95 was found to be 0.7 for regular patterns. Irregularities did not cause further deterioration. Five times layered repainting causes a statistically significant decrease in magnitude of the interplay effect across all breathing patterns by 50-80%, but is approximately 50% less effective against baseline shift than against other types of breathing. Interplay effect size correlates strongly with amplitude, but this correlation can be obscured because period and phase introduce very large variance. The interplay effect in general is large for investigated target size, prescription dose, beam configuration and machine performance. It can cause up to 100% of the CTV to receive a clinically unacceptable dose and lead to large inhomogeneities. Irregular breathing was not found to be notably worse. Repainting is effective, even against irregular breathing, but baseline shifts can undermine its effectiveness. Separately considering breathing irregularities for tumours similar to that investigated here is deemed of low importance; it is more important to properly model the magnitude of the interplay effect using an accurate, individualised breathing pattern.

Finding a favorable cure for cancer has been one of the main clinical challenges today. Nowadays the majority of the patients is treated with radiotherapy. Recently research in to a new paradigm of radiotherapy, so called FLASH radiotherapy, has opened up a new insight in to reducing negative side effects. FLASH dose rates (>40 Gy/s) have the ability to reduce the prevalence of negative side effects in patients without sacrificing tumor control. This study focuses on whether dose rate affects the degree of damage to healthy tissues caused by ionizing radiation. Understanding the biological mechanisms that influence the response to radiotherapy and specifically FLASH radiotherapy is therefore key knowledge in the development of new therapeutic protocols.This study evaluates dose-rate effects using zebrafish embryos as a small model organism. To study dose-rate effects the zebrafish are irradiated at 24 hours post fertilization and biomarkers are researched to asses dose-rate effects. They are irradiated using a Co60 irradiator providing a dose-rate of 0.015 Gy/s and a second one providing a dose-rate of 0.211 Gy/s. Several biomarkers for radiation damage are identified, including rate of embryonic development (hatching), DNA breaks and apoptosis. The DNA breaks are measured using a TUNEL assay and apoptosis through a caspase-3 assay. The results of the assay's are measured using confocal microscopy.Monte Carlo simulations of the experimental setup were performed to assess the uniformity of radiation and effective dose. They show a homogeneous dose throughout the sample, they show a 37% lower dose in water in the sample compared to air.The hatching biomarker shows that the hatching rates are 0.24 +-0.03 for non irradiated samples and 0.37 +- 0.31 and 0.43 +- 0.20 for samples that have been irradiated with respectively 0.008 Gy/s and 0.015 Gy/s with a total dose of 10 Gy. It shows that the hatching rates decreases for the irradiated embryos compared to the non irradiated samples however no significant difference can be found for different dose-rates. The sub-cellular biomarkers, the TUNEL assay and the caspase-3 assay, seem more promising in detecting dose-rate differences but more research is needed to find the dose-rate effects on DNA damage and apoptosis. In order to gain a fuller picture into the biological mechanisms of FLASH dose-rate radiation a greater variety of biomarkers and more research in to DNA damage and apoptotic biomarkers is needed.

Proton therapy efficiency can be described as the ratio between tumour and non-tumour dose, while the tumour receives the planned dose. This efficiency is limited by the energy deposition property of the proton. To enhance the efficiency beyond this physical limit, targeted nuclear reactions during proton therapy could be exploited. For this purpose the 11B(p,3a) reaction has been researched. The three alpha particles created have a high linear energy transfer (LET) and cause an increase in energy deposited at the reaction site. This reaction has a high cross section for low proton energies and protons have a low energy in the tumour area. Hence, proton 11B capture will more frequently occur in the tumour area, which can increase the tumour to non-tumour dose. Proton boron capture therapy (PBCT) has been studied, using Monte Carlo simulation and by conducting experiments, which suggested a 90 and 50% increase in energy deposited respectively. However, these results were contradicted by other research, stating the 11B(p,3a) reaction could not significantly increase energy deposited during proton therapy. This discussion formed the foundation of this thesis. Two Monte Carlo methods, MCNP6 and Geant4, were used to investigate the reproducibility of the results obtained in previous research. First, MCNP6 simulations were done, which did not result in significant dose increase for proton therapy. To investigate the reproducibility of MCNP6 a second Monte Carlo method, Geant4, was used. Geant4 simulations produced similar results to MCNP6, equally not reproducing the promising results of previous research. Neither Monte Carlo method included all possible reaction cross sections for alpha particle creation during PBCT. Therefore, a simple Boltzmann model was made to simulate proton and alpha particle transport, including these missing reaction cross sections. To simulate the effect of below 1 MeV proton boron capture, a non-dynamic Boltzmann model was used, which resulted in a non significant increase in alpha production during PBCT. Afterwards, an alpha-proton-alpha avalanche reaction was investigated. For the effects of this reaction to be investigated a dynamic Boltzmann model was used, which resulted in no significant increase in alpha production either. To investigate if the 11B(p,3a) reaction could cause the increase in cell death found by experimental work, the cell killing potential of each alpha particle created was calculated analytically. Using the alpha production rate from the non-dynamic model, this calculation resulted in ten cells killed per alpha particle created. The alpha particles created during PBCT have a track length similar to the size of one cell, therefore it is unlikely for an alpha particle to kill ten cell. None of the obtained results, provide evidence that the 11B(p,3a) reaction will enhance the efficiency of proton therapy. Thus, we hypothesize that the experimental results could be related to boron having a radio sensitizing effect on cells. To test this irradiating boronated cells with gamma rays is proposed as the next step for further research.

Hyperthermia is a form of thermal therapy used in combination with radiotherapy (RT) or chemotherapy (CTx) to enhance their effects. By using electromagnetic (EM) waves, tissue can be locally heated to temperatures between 39−44 ◦C. When applied to tumor tissue, this heating induces various biological responses that enhance both RT and CTx. The clinical workflow includes a hyperthermia treatment planning (HTP) stage, where treatment parameters are optimized to determine the phase and amplitude settings for the applicator antennas. HTP is subject to significant uncertainties. Variability in patient anatomy, dielectric and thermal tissue properties, patient positioning, and applicator modeling causes variations in the temperature distribution and thus affects the treatment. This study uses polynomial chaos expansion (PCE) to model these uncertainties and evaluate their impact on treatment parameters. PCE serves as a meta-model of the simulation software within a defined range of input variables. This meta-model allows for rapid calculation of temperature distributions, enabling the sampling of numerous scenarios to achieve statistical accuracy regarding the impact of treatment uncertainties. We successfully created a pipeline including patient modeling, treatment planning and uncertainty analysis using PCE. We developed separate smaller sub-models for different sources of uncertainty, including positioning, dielectric conductivity (σ), dielectric permittivity (ϵ), thermal conductivity (k), and perfusion rate (ω). Using these models we assessed the feasibility of building a PCE model for each uncertainty and determine the optimal settings for a comprehensive model combining input variables from different uncertainty sources. A high impact on the temperature achieved in 90% of the target volume (T90) was observed, with standard deviations observed up to 1.36 ◦C. Uncertainties in positioning, σ and ω had the largest impact, the impact of ϵ and k was significantly lower. The final model incorporated 31 parameters selected based on their impact on treatment parameters in a univariate analysis. This model was used to evaluate the combined effect of all treatment uncertainties, we did not succeed in maintaining the same accuracy as we show for the sub-models. We show that uncertainties are not necessarily additive, and that the combined effect of these uncertainties induced larger deviations in T90. The results demonstrated the potential of PCE to effectively handle complex uncertainty analysis in HTP, providing a robust framework for analyzing HTP. Future work should focus on validating these findings with a larger patient cohort to enhance the generalizability and reliability of the approach.

Combined chemo-radiation treatments are often used to enhance tumor cell kill, with respect to their separate treatments. Currently, chemo-radiation therapy treatments are prescribed as a separately optimized radiotherapy treatment with the addition of a generic chemotherapy regimen. The dosing of the combined therapies is mostly based on clinical experience. Hence there is a clear need for the optimization of dosing for combined chemo-radiation therapy treatments. The goal of this thesis is to propose novel approaches to optimize cancer treatment efficiency when combining chemotherapy and radiotherapy by using control-theoretical methods as well as radiation therapy models. A positive switched system is a dynamic system consisting of a set of Metzler state-space matrices and a specific switching law, which determines when and how to switch between the subsystems. Positive switched systems seem particularly appropriate to model cancer evolution under different treatments, and therefore to determine optimal treatment scheduling. They are used to model the effect of combined chemo-radiation treatments, thus enabling the systematic design of optimal therapy planning. The proposed model is especially appropriate for heterogeneous tumors, since it describes the impact of therapies on different tumor cell lines, and includes the possibility of mutations. Metastases are also included since the model considers separate tumor compartments, and the possibility of migrations between compartments are modeled as well. Undesired consequences due to toxicities are incorporated, and an upper bound can be set to limit the possible damage. A clinical dataset was obtained from the Erasmus MC, which was used to validate chemotherapy-alone model outcomes. Combined chemo-radiation treatments were also developed for a set of initial tumors. The model shows promising applications for cancer treatment design. The use of positive switched systems for combined chemo-radiation treatments is a new approach to obtain optimal drug usage and radiation fractionation schemes for each body compartment. The achieved findings could have a significant impact on effective treatment planning by choosing optimal drug durations or by using radiotherapy in combination with chemotherapy, for instance, to target certain chemotherapy-resistant areas, to achieve tumor control with minimal side effects.

In this report, a transformer based deep learning proton dose calculation algorithm called Dose Transformer Algorithm (DoTA) is described. This model learns to predict proton dose distributions by being trained with Monte Carlo generated data. Monte Carlo is the golden standard of proton dose calculation because it is very accurate, but it has relatively long computation times. In current proton therapy treatment programmes, Monte Carlo algorithms are the most commonly used models to perform dose calculation. The goal of the DoTA model is to predict proton dose distributions with Monte Carlo accuracy in the fraction of the computation time of Monte Carlo algorithms to speed up the dose calculation process in proton therapy treatment. The DoTA model can take patient geometry, random proton beam energy and random proton beam shape (2D Gaussian with different major and minor axes) as input. The addition of the random proton beam shape input is discussed in this report, together with a detailed explanation of the DoTA model. The DoTA model is trained with data from 9 lung cancer patients and 9 head & neck cancer patients. Being used on an Intel(R) Core(TM) i7-8565U CPU, the DoTA model managed to produce results with a gamma pass rate of 98.45 ± 2.60 % with an average computation time of 0.3 seconds. The gamma pass rate determines how similar the DoTA predicted dose is to the reference (Monte Carlo generated) dose. Compared to the average computation time of the Monte Carlo algorithm that was used to generate the training data, which is 20 seconds on the same CPU, we can conclude that the DoTA model has the potential to greatly improve dose calculation times in proton therapy treatment, especially when used on a system with greater processing power. Because the DoTA model is able to deliver accurate results in a small amount of time, it has the potential to be used for real-time dose calculation. Real-time dose calculation could account for small changes in patient geometry during treatment, which increases the accuracy of the treatment and minimizes side effects. The DoTA model can also be used for other radiotherapy types like phonon therapy and electron therapy (in that case it needs to be trained with different data).

In proton therapy, the calculation of the dose distribution is of great importance in order to find the best treatment plan. For a treatment plan with high quality, several error scenarios are investigated, in order to come up with a general plan that suits best in these scenarios. All of the scenarios have a different simulated error, e.g., in patient displacement and proton beam range. In order to cope with a large number of errors scenarios, a quick calculation method is needed. In this work, a method that used Reduced Order Modelling (ROM) and subsequently polynomial regression to overcome the computational problem is presented, as well as its results on a real head and neck cancer patient. Before error scenarios are simulated, using inverse optimisation a treatment plan is calculated for a perfectly known position and organ structure of the patient, i.e. the nominal scenario. Then, several error scenarios are simulated and for each of these scenarios, the real dose on tissue is calculated, using the treatment parameters that were previously set for the nominal scenario. The dose distributions of all these scenarios are stored in a matrix, on which a Singular Value Decomposition (SVD) is executed. A reduced order of singular values and vectors is used, together with regression models on the right singular vectors, to reconstruct the dose distribution matrix. This is done in order to create a computationally cheap method for dose distribution calculations. The right singular vectors are now a function of the used parameters: positioning errors dx, dy, dz and proton beam range error dρ. The dose distribution matrix is reconstructed from the left singular vectors, first few singular values and the function for the right singular vectors. Range errors with a uniform distribution with a maximum of 3% and positioning errors with a normal distribution with a standard deviation of 3 mm in each direction were simulated. For 100 error scenarios, the dose distribution matrix could be reconstructed with an average accepted error of 1% on a voxel relative to the maximum dose for 96:5% of the voxels, using a voxel-by-voxel comparison. These results were obtained for a 17th order of the SVD and a 7th order polynomial in the regression. On unseen test data, the acceptance was 81:0% of the voxels with an allowed error of 1%. For an allowed error of 3%, 92:5% of the voxels of this test set were accepted. These results are promising and encouraging for future research. First of all, it is recommended to examine the regression of the right singular vectors more closely, especially the regression of the higher order right singular vectors. Further subjects of interest would be to perform the same procedure that is done on the dose distribution matrix on the dose influence matrix. Also, harder-to-reach error scenarios, like organ movement and tumour deformation could be taken into account. This research was conducted at the Medical Physics and Technology section, a part of the Department of Radiation Science and Technology at Delft University of Technology. The project was part of the Bachelor program of Applied Physics at the Technical University of Delft.

In order to create radiotherapy treatment plans for cancer patients, dose calculations need to be done as quickly as possible to get accurate results. However, current dose calculation algorithms take too much time to be deployed effectively. The current in house algorithm of the Medical Physics and Technology Section at the TU Delft, attempts to solve this problem by utilising a deterministic algorithm that has a significant time advantage over Monte Carlo algorithms. However, this comes with the cost of inaccuracy, one of which is that it assumes all dose is deposited locally along the beam path. This is inaccurate as secondary particles created from non-elastic nuclear interactions can deposit their dose far from the beam path due to retaining significant kinetic energy. This thesis attempts to reduce this inaccuracy by mapping and quantifying the secondary particles to assess their contribution in non-local dose deposition. And analysing the relevant particle’s energy and angle distributions to gain insight into the development of the particle's characteristics with depth. Thereafter the relevant particle’s are then utilised as a source to emulate their production in a primary proton beam at different depths to obtain the relevant 3D dose distributions. The analysis concluded that secondary protons are the most relevant secondary particle as they contribute to 88% of the secondary dose and have a significant range to deposit their dose non locally. By utilising the secondary protons as a source, it was found that the relative error between the integrated depth dose (IDD) of the scored protons and the IDD obtained directly from Monte Carlo simulations is equal to 5.1% in the z-direction and 3.4% in the x and y-direction. The absolute difference was found to be 1.54 × 10−5 Gy which is equal to 0.096% of the total dose and 2.75% of the dose contributed by all secondary particles. The results show that the methodology can produce accurate 3D dose matrices for secondary protons at different depths, which can then be used to improve the accuracy of the in house algorithm by adding the precalculated 3D dose matrices to the algorithm

Numerical models are paramount in describing the complex physical world around us. They are often based on non-linear functions, which have to be solved in order to run a simulation. The first goal of this thesis is to compare the Jacobian-free Newton-Krylov method against the regular Newton-Raphson method for solving non-linear equations. When doing this, preconditioning was only used for the Newton-Raphson method. As a second goal, the optimal parameters for solving a coupled system with the Schur complement method are determined. For computational purposes, these analyses are performed using PETSc. The comparison between JFNK and Newton-Raphson were performed by solving a heat equation. The physical system which was analyzed is a one-dimensional radiating rod, with a heterogeneous thermal conductivity and Dirichlet boundary conditions. Firstly this rod was modelled as one single system. For this case, it was found that Newton-Raphson outperformed the JFNK method by a factor of 5-1700, depending on the type of preconditioning used. Secondly, the rod was modelled as a coupled system by solving two parts of the rod separately. In this case, it was found that the JFNK method outperformed the preconditioned Newton-Raphson method by a factor of 9.1± 0.3. It was concluded that the JFNK is only favorable over regular Newton-Raphson for coupled systems. To achieve the second goal, an incompressible Navier-Stokes coupled system was solved using the Schur complement method. The coupled system originated from incompressible flow in a back-step pipe, using a finite element method. For the Schur complement parameters, it was shown that using an approximation of the Schur complement offered a significant increase in efficiency. The momentum and pressure subsystems were analyzed separately. The momentum subsystem performed best with a relative tolerance of τr =10-4.5, while the pressure subsystem performed best with a tolerance of 10-2. With these parameters, the Schur complement method had the same computational time as the pressure-correction method, which was used as a benchmark. For future research, there are three main recommendations. Firstly, to use preconditioning for the JFNK method. This was not done in this thesis because of technical constraints, but it is theoretically possible. Secondly, if the JFNK method is preconditioned, it could be used for the incompressible Navier-Stokes simulation. Finally, there are combinations of parameters which were not tested for the Schur complement method. Trying out more combinations might result in finding an even more optimized method.