finding norm of vector but vector not visible in workspace

function pdex1
clear;
%PDEX1  Example 1 for PDEPE
%   This is a simple example with an analytical solution that illustrates
%   the straightforward formulation, computation, and plotting of the solution
%   of a single PDE.
%
%   In the form expected by PDEPE, the single PDE is
%
%   [1] .*  D_ [u] =  D_ [   Du/Dx   ] +      [0]
%              Dt        Dx
%    ----         ---       -------------    -------------
%     c            u        f(x,t,u,Du/Dx)   s(x,t,u,Du/Dx)
%
%   The equation is to hold on an interval 0 <= x <= 1 for times t >= 0.
%
%   Two kinds of boundary conditions are chosen so as to show how they
%   appear in the form expected by the solver.  The left bc is u(0,t) = 0:
%
%     [u]    +     [0] .* [       Du/Dx          ] = [0]
%
%     ---          ---    ------------------------   ---
%   p(0,t,u)      q(0,t)        f(0,t,u,Du/Dx)        0
%
%   The right bc is taken to be u(1,t) + Du/Dx(1,t) = 0:
%
%   u(1,t) + [1] .* [   Du/Dx    ] = [0]
%
%   --------------   ---    --------------   ---
%      p(1,t,u)     q(1,t)   f(1,t,u,Du/Dx)   0
%
%   The problem is coded in subfunctions PDEX1PDE, PDEX1IC, and PDEX1BC.
%
%   See also PDEPE, FUNCTION_HANDLE.

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2014 The MathWorks, Inc.

%we are solving for u where u is T - T_infinity (T is temperature)at any
%point x along the length of the cantilever

global rho Cp K h length w d A sigma beam_size P0 B D x0 total_time time_steps length_steps omega time_at_eval p q position_eval lambda
length = 200;  %cantilever length, microns
w = 20; % width of cantilever, microns
d = 5; %thickness of cantilever, microns
A = w*d;  %cross sectional area of cantilever, micron^2
x0 = (1/4)*length;  %illumination spot location as a fraction of cantilever length, microns
h = 10*10^-12;  %convection coefficient, W/um^2-K
K = 148*10^-6;  %thermal conductivity of cantilever material, W/um-K
rho = 2.33*10^-15;  %density of cantilever material, kg/(um)^3
Cp = 700; %specific heat of cantilever material, J/kg-K
beam_size = 5; %laser spot size, microns
sigma = beam_size/6; %standard deviation as a function of laser spot size, microns
P0 = 0.5*10^-3; %applied input power, Watts
omega = [10000 50000]; %laser modulation frequency, kHz
B = P0/(((2*pi)^0.5)*sigma*K*A);  %coefficient of gaussian term, kelvin/(um)^2
D = rho*Cp/K; %coefficient of transient term, second/(um)^2
total_time = 0.01;  %total time duration for the simulation
time_steps = 500; %time steps taken
length_steps = 500;  %number of steps in length x
p = 1/4; %fraction of total duration time; used to set intermediate time at which we want to evaluate spatial temperature distribution
q = 1/4; %fraction of total length; used to set any intermediate length value at which we want to evaluate temporal temperature distribution
time_at_eval = total_time*p; % time at which we want to evaluate spatial temperature distribution
position_eval = length*q;
lambda = K*time_steps/(length_steps);
m = 0; %for cartesian coordinates
x = linspace(0,length,length_steps);  %defining x mesh vector along length of cantilever
t = linspace(0,total_time,time_steps);  %defining time mesh vector
disp(omega);

for k=1:2
    pdex1 = @(x,t,u,DuDx) pdepar(x,t,u,DuDx,omega(k));
    sol{k} = pdepe(m,pdex1,@pdex1ic,@pdex1bc,x,t);
    u = sol{k}(:,:,1);
    a = norm(u);

% A simple X-Y plot of temperature as a function of distance is generated below at a particular time
figure(5);
hold on
plot(x,u(round(end*p),:),'-o');
title(['temperature profile along x at t = ', num2str(time_at_eval), ' sec']);
legend(strcat(num2str(omega(1)), ' kHz'), strcat(num2str(omega(2)),' kHz'),'Location', 'Northeast');
xlabel('Distance x (um)');
ylabel('temperature difference wrt surrounding (deg K)');

% A simple X-Y plot of temperature as a function of distance is generated below at a particular time
figure(6);
hold on
plot(t,u(:,position_eval),'-o');
title(['Temperature profile vs time at x = ', num2str(position_eval), ' um']);
legend(strcat(num2str(omega(1)), ' kHz'), strcat(num2str(omega(2)),' kHz'), 'Location', 'NorthEast');
xlabel('Time t (sec)');
ylabel('temperature difference wrt surrounding (deg K)');
end
% pdex1 = @(x,t,u,DuDx) pdepar(x,t,u,DuDx, omega);
% sol = pdepe(m,pdex1,@pdex1ic,@pdex1bc,x,t);
% Extract the first solution component as u.  This is not necessary
% for a single equation, but makes a point about the form of the output.


% A surface plot is generated below.
% figure;
% surf(x,t,u);
% title(['Numerical solution computed with ', num2str(length_steps), ' mesh points.']);
% xlabel('Distance x');
% ylabel('Time t');


% A simple X-Y plot of temperature as a function of distance is generated below at time t = 2 secs.
% figure;
% plot(x,u(end,:),'-o');
% title(['temperature profile along x at t = ', num2str(total_time), ' sec']);
% legend('Numerical','Analytical', 'Location', 'NorthEast');
% xlabel('Distance x (um)');
% ylabel('temperature difference wrt surrounding(deg K)');

% A simple X-Y plot of temperature as a function of distance is generated below at time t = 2 secs.
% figure;
% plot(t,u(:,x0),'-o');
% title(['Temperature profile vs time at x = ', num2str(x0), ' um']);
% legend('Numerical','Analytical', 'Location', 'NorthEast');
% xlabel('Time t (sec)');
% ylabel('temperature difference wrt surrounding (deg K)');
% --------------------------------------------------------------------------

% function [c,f,s] = pdex1pde(x,t,u,DuDx)  %main PDE solver
% global B D x0 sigma omega k
% c = D;   %coefficient of dT/dt
% f = DuDx;   %at f = DuDx, we get double derivative d2u/dx2
% for k=1:7
%     s(k) = B*(1+0.75*sin(omega(k)*t))*exp(-(x-x0)^2/2*(sigma)^2)   %gaussian source term with sinusoidally varying power
% end

function [c,f,s] = pdepar(x,t,u,DuDx,omega)
global B D x0 sigma
c = D;   %coefficient of dT/dt
f = DuDx;   %at f = DuDx, we get double derivative d2u/dx2
s = B*(1+0.75*sin(omega*t))*exp(-(x-x0)^2/2*(sigma)^2);   %gaussian source term with sinusoidally varying power


% --------------------------------------------------------------------------


function u0 = pdex1ic(x)  %initial condition function
u0 = 0;  %initial temperature difference is zero since the cantilever is at room temperature at t = 0

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)   %boundary condition function. every boundary condition needs to be of the form p + q*f. L and R are subscripts for left and right face
global h K
pl = ul; %temperature difference at left end or clamped end of cantilever is zero since it is kept at room temperature
ql = 0;
pr = ur*(h/K); %convective boundary condition at right end or free end of cantilever
qr = 1;

[Akshay Deolia]

Forgot to add..I'm trying to get the norm of the solution vector for the time-varying temperature solution. The time-varying temperature solution is not converging and hence I'm trying to find the L2 norm of the solution vector for a given number of time mesh points.

[Stephen23]

What is the point of clearing an empty function workspace?

Answer 1

[Sulaymon Eshkabilov]

You'd need to assign an output or outputs of variables computed within your fcn file, e.g:

function u = pdex1
...
    

That delivers you all computed u data into workspace.

[Akshay Deolia]

Hi Sulaymon,

The solution is stored in a 50x100 matrix for 50 time steps and 100 length steps. Is is possible for me to store this variable? For example, assigning it to a variable "a" and then I'd like to extract a particular column or row after that.

[Sulaymon Eshkabilov]

yes, of course. It is a very simple:

function [u, a] = pdex1
...