Home /
Expert Answers /
Computer Science /
below-is-my-code-error-that-i-see-while-running-the-tutotial-main-m-file-error-using-diff-dif-pa617

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<maxiter) % iterate until one of the criteria is met

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 < desired tolerance

elseif norm(f) <= tol

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=?_{r_{i}}r1?(?_{??}??_{rr})dr$ where $P_{i}$ is the luminal pressure applied to the vessel (assuming $P_{o}?0$ ), $r_{i}$ and $r_{o}$ 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 $R_{i}$, the outer radius in the reference configuration $R_{o}$, the axial stretch $?_{z}$ in the current configuration, the luminal pressure $P_{i}$ in the current configuration, and the constitutive parameters. You will provide an initial estimate of the outer radius in the current configuration $r_{o}$ 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=r_{i}+?_{z}1?(R_{2}?R_{i})?$ 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=?_{r_{i}}r1?(?_{??}??_{rr})dr$ as the right hand side of Equation 1. Your goal is to find the value of $r_{o}$ that minimizes the difference between $T$ and $P_{i},f=T?P_{i}$, 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=?R?r??Rr???_{s}=1$, from which: $dr=?_{z}?rR?dR.$ 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 $P_{i}$ 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 $?_{r_{1}}r1?(?_{n^?^}??_{rr})dr$ through the thickness, from the inner radius $r_{i}$ (given) to the outer 1. Set a toleranoe and nuximun number of iterations to represent your criteria to end the root radius $r_{o}$ (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 $R_{i}$ to $R_{o}$. 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 $(x_{0})$ and evaluate $f$ using the MATI.AB function fecoll(). 1. Define a vector of radii $R?[R_{i}:$ step $:R_{u}]$ in the refcrence configuration. Choose the step 4. To define the incremental step, set $xinc$ to $x_{n}+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 $?_{r_{1}}r1?(?_{00}??_{r7})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 $R_{1}$ and $R_{2}$ to use for the call call the function $f$ for both steps. numerical integration. 6. Iterate with the Newton-Raphson, $x=x_{0}?ff?$, method to find the next $x$. (b) Calculate $r_{1}$ and $r_{2}$ from $R_{1}$ and $R_{2}$, 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 $F_{1}$ and $F_{2}$ at radii $r_{1}$ and $r_{2}$, 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 $(r_{0})$ as the root of equation 1. $????r=r(R)?=?z=?_{z}Z?$ Express $F$ as a $3×3$ 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 $3×3$ diagonal matrix. The strain energy function for the Nec-Hookean model is defined as $W?c_{1}(I_{1}?3)$ which $c_{1}$ is a material purameter and $I_{1}$ 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 $_{r}(?_{?g}??_{rr})dr$ at $R_{1}$ and $R_{2}$ 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 $R_{1}$ and $R_{2}$, you will also nexd to express

INTRODUCTIONThe Newton-Raphson approach starts with a pre