Below is my code

Error that i see while running the Tutotial_main.m file

Error using diff Difference order N must be a positive integer scalar.

Error in Newton_Raphson_tutorial (line 35) f_prime0 = diff(f,x0,xinc); % compute the derivative of f, between x0 and xinc

Error in Tutorial_m (line 51) x = Newton_Raphson_tutorial(H,x0); % call the Newton Raphson function (Newton_Raphson_tutorial.m)

for Tutorial_main.m

%=========================================================================

% Lecture 16: In Class Tutorial

%

% This function calculates the radial equilibrium function for an axially

% stretched and pressurized thick wall vessel and is part of the set of

% equations you will implement for your vasculature project

%

% Input data:

% luminal pressure (Pi), axial stretch (lambdaz_v)

% material parameters, radii in ktf (Ri, Ro)

%

% Output data:

% approximation of the outer radius, ro

%

% The inverse solution of the radial equilibrium involves finding

% the root of the equation:

% Pi - int_{ri}^{ro} (tqq-trr)/r dr = 0

%=========================================================================

%cd YOUR DIRECTORY;

close all, clear all, clc

%% Load and Define INPUT DATA

global input_data

% Conversion factor for Pressure: mmHg to kPa

% mmHg2kPa = 101.325/760;

% Load results from NLR

input_data = load('output_NLR_WT_ATA_FIT_obj5_a1.mat');

% Pressure & Axial Stretch @ loaded configurations of interest

lambdaz = input_data.data_kl.lambdaz_exp; %stored in input data

Pi = input_data.data_kl.Psys_exp; %in-vivo pressure

%% Solve the equilibrium equation

% Equilibrium equation for a pressurized and axially stretched vessel

% Unknown: outer radius, ro

% ro will be estimated as the roots of the equilibrium equation

H=cell(1,1);

H{1}=@equilibrium_r_or_loaded_tutorial; %handle to equilibrium equation

% Parameter to be estimated in loaded configuration is the outer radius (ro)

% x0 is the initial guess for ro based on the input experimental data

x0 = [input_data.data_ktf.or_exp];

input_data.data_kl.lambdaz = lambdaz; % axial stretch

input_data.data_kl.Pi = Pi; % liminal pressure, in kPa

% Newton Raphson Scheme, Inputs: function handles and evaluation points

x = Newton_Raphson_tutorial(H,x0); % call the Newton Raphson function (Newton_Raphson_tutorial.m)

or_est = x(1,1); % calculated outer radius

CODE for Newton_Raphson_tutorial.m

Doubts

According to the instructions, I have tried to define the variables.

I am not sure what should be the tolerance(tol ), step size(h).

How to define f_prime0.

2.3 point 3 to evaluate f using feval() I don't know where to write that and how to write that.

function [x] = Newton_Raphson_tutorial(H,x0)

%=========================================================================

% Newton Raphson Scheme, X_n = X_(n-1) - f/f'

%

% H is the handle to the equilibrium equation

% x0 is the initial guess for the outer radius

%

% x is the approximate solution of the Newton Raphson method (in this case,

% ro)

%=========================================================================

%% Iteration Constraints

tol = 2; % tolerance criteria to stop scheme

maxiter = 4; % maximum number of iterations allowed before stopping

%% Newton Raphson Iteration

if length(H) ~= length(x0) % check if the number of equations match the number of unknowns

disp('1Error: number of unknowns does not coincide with number of equations')

return

else

% first perform Newton Raphson for the initial guess

% evaluate the function, f, with the initial guess, x0

f = equilibrium_r_or_loaded_tutorial(x0); % call: equilibrium_r_or_loaded_tutorial

iter = 0; % initial iteration

h = 0.01; % step size to perform the derivative (forward difference)

xinc = x0; % initial guess of the root

xinc = xinc + h; % current root + derivative step

f_prime0 = diff(f,x0,xinc); % compute the derivative of f, between x0 and xinc

u0 = f/f_prime0; % divide f by its derivative, f/f'

x = x0-u0; % update guess for ro: X = X0 - f/f' (Newton-Raphson Method)

% Now, continue iterating until either a solution is found (within the

% defined tolerance) or you have reached the maximum number of

% iterations

while (iter

iter = iter+1; % current iteration number

f = feval(f,x); % evaluate f at x

xinc = x; % current root guess

xinc = x + h; % current root + derivative step

f_prime = diff(f,x,xinc); % compute the derivative of f, between x and xinc

u = f/f_prime; % divide f by its derivative, f/f'

x = x0 - u; % update guess for ro: X_n = X_(n-1) - f/f' (Newton-Raphson Method)

end

end

%% Outputs

% Stop when we meet the max number of iterations

if iter == maxiter

disp ('---')

disp ('Maximum number of iterations reached')

disp ('---')

% or stop when solution is found: the error

elseif norm(f)

disp ('---')

disp ('Minimum tolereance reached')

disp ('---')

end

iter % number of iterations you performed

f % function evaluated

x % root, ro

end

CODE FOR equilibrium_r_or_loaded_tutorial.m

function [fvalue] = equilibrium_r_or_loaded_tutorial(x0)

%=========================================================================

% RADIAL EQUILIBRIUM - Local

%

% This function calculates:

% Pi - int_{ri}^{ro} (tqq-trr)/r dr = 0

%

% x0 is the current solution of the Newton-Raphson method

%=========================================================================

global input_data

ro = x0(1,1); % outer radius

Ri = input_data.data_ktf.ir_exp; % reference inner radius

Ro = input_data.data_ktf.or_exp; % reference outer radius

lambda = input_data.data_kl.lambdaz; % axial stretch

Pi = input_data.data_kl.Pi; % in vivo pressure

ri = sqrt(ro.^2-1./lambda*(Ro^2-Ri^2)); % calculate the inner radius

a = Ri; % lower limit of the independent variable a

b = Ro; % upper limit of the independent variable b

T = 0; % initial value for the integral (T is the result of the integral)

n = 10; % number of spatial steps

h = 0.1; % spatial step size, based on n and the bounds [a,b]

for index1 = 0:n-1 % for loop going from inner to outer radius

R1 = Ri+index1*h; % reference position

r1 = sqrt(ri.^2+1./lambda*(R1^2-ri^2)); % calculate the radius by mapping from the reference configuration

F1 = diag([ r1 r1 lambda ]); % define the the deformation gradient tensor (use the diag function)

sigma_extra1 = constitutive_model_NH(F1); % calculate sigma_extra using the constitutive model

f1 = (R1/lambda*r1); % evaluate the function inside the integral, at R1

R2 = R1+h; % next position

r2 = sqrt(ri.^2+1./lambda*(R2^2-ri^2)); % calculate the radius by mapping from the reference configuration

F2 = diag([ r1 r1 lambda ]); % define the the deformation gradient tensor (use the diag function)

sigma_extra2 = constitutive_model_NH(F2); % calculate sigma_extra using the constitutive model

f2 = (R1/lambda*r1); % evaluate the function inside the integral, at R2

% calculate the intergral by solving the Trapezoidal Rule

T = f1+f2; % remember to add the previous T_(index-1)

end

fvalue = T - Pi; % ultimately fvalue needs to go to zero (within the defined tolerance)

end

Use Numerical Differentiation, Integration, and Root Finding to Solve Radial Equilibrium for a Blood Vessel Learning Objectives By completing this in-class MATLAB tutorial you will learn how to: 1. resolve radial equilibrium for an axially-extended and pressurized thick-wall vessel; 2. implement the Newton-Raphson algorithm for root finding; 3. implement a numerical differentiation scheme; 4. implement a numerical integration scheme. 1 Radial Equilibrium The equation you will be solving in this tutorial (and for your vascular project!) is the radial equilibrium for an axially-extended and pressurized thick wall vessel: Pi=riror1(rr)dr where Pi is the luminal pressure applied to the vessel (assuming Po0 ), ri and ro are the inner and outer radii in the current configuration, respectively, and rr(r) and (r) are the components of the Cauchy stress in the radial and circumferential direction, respectively. You are given input data, which includes the inner radius in the reference configuration Ri, the outer radius in the reference configuration Ro, the axial stretch z in the current configuration, the luminal pressure Pi in the current configuration, and the constitutive parameters. You will provide an initial estimate of the outer radius in the current configuration ro and you will refine the estimate iteratively, by implementing the Newton-Raphson method. You will predict the radii in the current configuration using equation 2, which maps each radius R in the reference configuration into the radius r in the current configuration. r=ri2+z1(R2Ri2) You will predict the component of the Cauchy stress at each material point through the thickness using a given constitutive model. 2 Root Finding Let's define T=riror1(rr)dr as the right hand side of Equation 1. Your goal is to find the value of ro that minimizes the difference between T and Pi,f=TPi, by using the Newton-Raphson method (see lecture 15), which requires you to evaluate f as well as its derivative, f. To evaluate f you will need to use numerical integration. Page 1 2.1 Input Data Assignment 1. Open Tutorial_main.m, load the Matlab structure output_NLR_WT_ATA_FIT_obj5_a1.mat, and arsign it to the global variahle moutanta. dr as a function of dR. To do so, recall that detF=RrRrs=1, from which: dr=zrRdR. 2. The axial stretah is stored into variable intud data.daha blandwhaz exp. Assign it to variable lambdaz. (f) Sum the contributions at all iterations to obtain an estimate of the integral T. 3. The luminal pressure (in kPa ) is stored into variable inpucl-dala_dufa-dl.Psys_exp. Assign to 2.3 Newton-Raphson method vuriable Pi You will approximate the outer radius with the Newton Raphson method. To do so, first open 2.2 Numerical Integration Nentor_Raphom_tutorial.m. Open equibibiun r or loweled futorial.m. You will obtain T frun Feguation 3 by numerically integrating r10r1(n^^rr)dr through the thickness, from the inner radius ri (given) to the outer 1. Set a toleranoe and nuximun number of iterations to represent your criteria to end the root radius ro (guessexd, to be optimised). To do s), you will nexd to calculate the radius in the current finding process. contiguration r for each value of the reference radius R, from Ri to Ro. Furthermore, you will need 2. Set your step size for the derivative (h). Lagrange multipliur (sigmane-rfra) at mach radial position r. 3. Start, with your initial guess (x0) and evaluate f using the MATI.AB function fecoll(). 1. Define a vector of radii R[Ri: step :Ru] in the refcrence configuration. Choose the step 4. To define the incremental step, set xinc to xn+h (for the initial iteration) and xinc to x+h size to ensure accuracy of the numerical integration algorithm. (for the rest of the iterations). This will be used to find the derivative of f. 2. Use a for loop to numerically integrate function r1r0r1(00r7)dr using the trapezoid rule. 5. Fvaluate the derivative my imphementing a forwanl difference numerical sheme. Note, you (a) for each iteration of the for laop identify two adjacent radii R1 and R2 to use for the call call the function f for both steps. numerical integration. 6. Iterate with the Newton-Raphson, x=x0ff, method to find the next x. (b) Calculate r1 and r2 from R1 and R2, respectively, using Equation 2 . 7. Continue with the while function to itarate on Newton-Raphson methoxd until either criteria (c) Calculate the deformation gradient tensors F1 and F2 at radii r1 and r2, respectively, is met (satisfied the tolerance criteria or reached maximum number of iterations). using the mapping for an axially-cxtendexd and pressurized vesel. Call the Newton_Raphison_tutorial function in Tutorial_main.m to calculate outer radius (r0) as the root of equation 1. r=r(R)=z=zZ Express F as a 33 diagonal matrix containing the non-zero, diagonal components of the deformation gradient tensor. (d) Calculate the extra, deformation-dependent components of the Cauchy stress tensor using a Neo-Ilookean constitutive model and assign them to variable signa_ertra, which is also a 33 diagonal matrix. The strain energy function for the Nec-Hookean model is defined as Wc1(I13) which c1 is a material purameter and I1 is the first invariant of the right Cauchy-Green deformation tensor C. To calculate the stress, use the provided Mat lab function censtifution-medel.NH.m, which receives the deformation gradicnt tensor F as an input and provides the signa-extra as output. (e) Evaluate r1(grr)dr at R1 and R2 and calculate their contribution to the throughthickness integral. Because you are expressing current radii r and Cauchy stress compoments sigma -erlra using the reference variables R1 and R2, you will also nexd to express