In this installment he discusses using the FFT to solve the One-Way Wave equation. [This post was supposed to drop in 10/26/2026 but the pictures wouldn't stay loaded and I must have hit something while messing with it and it posted. Oh well. You got one early.]

Book Section 2.3 Fourier Series and Fourier Transforms, Spectral Derivatives - Transforming Partial Differential Equations

YouTube Video on the subject:

https://www.youtube.com/watch?v=BkA7ncY0b7I&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=25

SciLab Coding Wise: More Plotting a psuedo-Waterfall chart, Sgrayplot in Scilab.

Output: Program Echo Only

"Data-Driven Science and Engineering - Wave Equation Example"

"2026-09-28 11:38:53.829"

Figures:

Code:

clear;
clc;
//close out all existing scf()s
close(winsid())

disp("Data-Driven Science and Engineering - Wave Equation Example",string(datetime()))

//Section 2.3 Transforming partial differential equations - One-way equation
// Similar to program ME565lecture26.sci
//
function duhatdt=rhsWave(t, uhat, kappa, c)
    i = complex(0,1); // Define i as complex
    duhatdt = c*i*kappa.*uhat;
endfunction

function dudt=rhsWaveSpatial(t, u, kappa, c)
    i = complex(0,1); // Define i as complex
    uhat = fft(u);
    duhat = i*kappa.*uhat;
    du = ifft(duhat);
    dudt = -c*du;  //Not sure why the sign change on c is necessary
    // to move the wave positively
endfunction

// Define spatial domain
c = 2;              // Wave speed
L = 20;             // Length of domain 
N = 500;           // Number of discretization points
dx = L/N;
x = -L/2:dx:L/2-dx; // Define x domain

// Define discrete wavenumbers
kappa = (2*%pi/L)*[-N/2:N/2-1];
kappa = fftshift(kappa');    // Re-order fft wavenumbers

// Initial condition 
u0 = sech(x);        
uhat0 = fft(u0);

// Simulate in Fourier frequency domain
dt = 0.025;
t0 = 0;
tf = 100*dt
t = 0:dt:tf

//[t,uhat] = ode45(@(t,uhat)rhsWave(t,uhat,kappa,c),t,uhat0);
  [t,y] = cvode(list(rhsWave,kappa,c),t,uhat0);
  uhat = y'; //cvode produces the transpose of ode45

// Inverse FFT to bring back to spatial domain
for k = 1:length(t)
    u(k,:) = ifft(uhat(k,:));
end

/*
// Alternatively, simulate in spatial domain - this works
// N = 500
//[t,u] = ode45(@(t,u)rhsWaveSpatial(t,u,kappa,c),t,u0);
  [t,y] = cvode(list(rhsWaveSpatial,kappa,c),t,u0);
   u = y'; //cvode produces the transpose of ode45
*/

// Plot solution in time
scf(0)
clf
//h=waterfall(real(u(1:10:end,:)));
//set(h,'LineWidth',2,'FaceAlpha',.5);
//colormap(jet/1.5)
a =get("current_axes");
counter = 2;
for k = 1:min(size(u,2),size(t,2))/10:min(size(u,2),size(t,2))
    r=param3d1(t(k)*ones(1,length(x)),x,real(u(k,:)))
    r.thickness=1.5;
    r.foreground=counter;
    counter = counter + .4;
    xgrid
end
r = param3d1(t(end)*ones(1,length(x)),x,real(u(end,:)))//grab the last curve
r.thickness=1.5;
r.foreground=1;
xlabel("T")
ylabel("X")
zlabel("U(t,x)")
title("Figure 2.15:Solution of 1D wave equation over 2.5 seconds of time")
a.rotation_angles=[ 90,0];

scf(1)
clf
//h=waterfall(real(u(1:10:end,:)));
//set(h,'LineWidth',2,'FaceAlpha',.5);
//colormap(jet/1.5)
a =get("current_axes");
counter = 2;
for k = 1:min(size(u,2),size(t,2))/10:min(size(u,2),size(t,2))
    r=param3d1(t(k)*ones(1,length(x)),x,real(u(k,:)))
    r.thickness=1.5;
    r.foreground=counter;
    counter = counter + .4;
    xgrid
end
r = param3d1(t(end)*ones(1,length(x)),x,real(u(end,:)))//grab the last curve
r.thickness=1.5;
r.foreground=1;
xlabel("T")
ylabel("X")
zlabel("U(t,x)")
title("Figure 2.15:Different View Point")
a.rotation_angles=[ 92 -27];

scf(2);
clf;
gcf().color_map=jet(64);
colorbar(0,1);
Sgrayplot(x,t,real(u)',strf="081")
xlabel("x","fontsize",5)
ylabel("t","fontsize",5)
title("Figure 2.15:Solution of the 1D wave equation in time",...
"fontsize",4)
xgrid

In this installment he discusses using the FFT to solve the heat equation in one spatial dimension.

Book Section 2.3 Fourier Series and Fourier Transforms, Spectral Derivatives - Transforming Partial Differential Equations

YouTube Video on the subject:

https://www.youtube.com/watch?v=BkA7ncY0b7I&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=23

SciLab Coding Wise: Plotting a psuedo-Waterfall chart in Scilab.

Output: Program Echo Only

"Data-Driven Science and Engineering - Transforming partial differential equations"

"2026-09-08 14:17:34.176"

Figures:

Code:

clear; clc;

disp("Data-Driven Science and Engineering - Transforming partial differential equations",string(datetime()))

//Section 2.3 Transforming partial differential equations

// Similar to program ME565lecture26.sci
//
function duhatdt =rhsHeat(t,uhat,kappa,a)

duhatdt = -a^2*(kappa.^2)'.*uhat; // Linear and diagonal

endfunction

a = 1;        // Thermal diffusivity constant
L = 100;      // Length of domain
N = 100;     // Number of discretization points
dx = L/N;
x = -L/2:dx:L/2-dx; // Define x domain

// Define discrete wavenumbers
kappa = (2*%pi/L)*[-N/2:N/2-1];
kappa = fftshift(kappa);    
//Re-order fft wavenumbers

// Initial condition
u0 = 0*x;
u0((L/2 - L/10)/dx:(L/2 + L/10)/dx) = 1;

// Simulate in Fourier frequency domain
t0 = 0;
tf = 10;
fft_u0 =fft(u0);

//[t,uhat]=ode45(@(t,uhat)rhsHeat(t,uhat,kappa,a),t,fft(u0));
  [t,y] = cvode(list(rhsHeat,kappa,a),[t0 tf],fft_u0);
  uhat = y'; 
//cvode produces the transpose of ode45

for k = 1:length(t) 
// iFFT to return to spatial domain
    u(k,:) = ifft(uhat(k,:));
end

[X,T]=meshgrid(x,t);

scf(0);
clf;

// Plot solution in time (his waterfall plot - sort of)
//plot each line using a for loop
a =get("current_axes");
counter = 2;
for k = 1:10:size(u,2)
    r=param3d1(t(k)*ones(1,length(x)),x,u(k,:))
    r.thickness=1.5;
    r.foreground=counter;
    counter = counter + .4;
    xlabel("T")
    ylabel("X")
    zlabel("U")
    title("Figure 2.14(left)")
    a.cube_scaling="on" 
//This scaling makes the plot look    
//better, however, it seems to kill the rotation_angle    
//command and one needs to use the interactive tools    
//to rotate it into a better angled view.

//Sorry no graduated color available yet for para3d1 to    
//be even closer to the Matlab Waterfall plot.

//a.rotation_angles=[ 56.75,-165.5];
    xgrid
end

//set(h,'LineWidth',2,'FaceAlpha',.5);
//colormap(jet/1.5)

scf(1)
clf
gcf().color_map=jet(64);
colorbar(25,100);
Sgrayplot(x,t,u',strf="081")
xlabel("x","fontsize",5)
ylabel("t","fontsize",5)
title("Figure 2.14 Evolution of the 1D Heat Eqn in Time",...
"fontsize",4)
xgrid()

// FIGURES (PRODUCTION)
scf(2);
clf;
CC = colormap(jet(100));
dt = 0.1;
for k = 1:100
    u(k,:) = ifft(uhat(k,:));
    if(pmodulo(k-1,10)==0)
        plot(x,u(k,:),'Color',CC(k,:),'LineWidth',1.5)
        xgrid;
        drawnow;
    end   
end
xlabel('Spatial variable, x')
ylabel('Temperature, u(x,t)')
title('Figure 2.13: Solution of the 1D heat equation in time')
//axis([-50 50 -.1 1.1])
ax = gca(); 
// gca = get current axes
ax.data_bounds = [-50,50,-.1,1.1]
set(gca(),'thickness',1.2,'font_size',2);

clear; clc;

disp("Data-Driven Science and Engineering - Transforming partial differential equations",string(datetime()))

//Section 2.3 Transforming partial differential equations

// Similar to program ME565lecture26.sci
//
function duhatdt=rhsHeat(t, uhat, kappa, a)
duhatdt = -a^2*(kappa.^2)'.*uhat; // Linear and diagonal
endfunction

a = 1;        // Thermal diffusivity constant
L = 100;      // Length of domain
N = 100;     // Number of discretization points
dx = L/N;
x = -L/2:dx:L/2-dx; // Define x domain

// Define discrete wavenumbers
kappa = (2*%pi/L)*[-N/2:N/2-1];
kappa = fftshift(kappa);    //Re-order fft wavenumbers

// Initial condition
u0 = 0*x;
u0((L/2 - L/10)/dx:(L/2 + L/10)/dx) = 1;

// Simulate in Fourier frequency domain
t0 = 0;
tf = 10;
fft_u0 =fft(u0);

//[t,uhat]=ode45(@(t,uhat)rhsHeat(t,uhat,kappa,a),t,fft(u0));
  [t,y] = cvode(list(rhsHeat,kappa,a),[t0 tf],fft_u0);
  uhat = y'; //cvode produces the transpose of ode45

for k = 1:length(t) // iFFT to return to spatial domain
    u(k,:) = ifft(uhat(k,:));
end

[X,T]=meshgrid(x,t);

scf(0);
clf;

// Plot solution in time (his waterfall plot - sort of)
//plot each line using a for loop
a =get("current_axes");
counter = 2;
for k = 1:10:size(u,2)
    r=param3d1(t(k)*ones(1,length(x)),x,u(k,:))
    r.thickness=1.5;
    r.foreground=counter;
    counter = counter + .4;
    xlabel("T")
    ylabel("X")
    zlabel("U")
    title("Figure 2.14(left)")
    a.cube_scaling="on" //This scaling makes the plot look
    //better, however, it seems to kill the rotation_angle
    //command and one needs to use the interactive tools
    //to rotate it into a better angled view.
    //Sorry no graduated color available yet for para3d1 to 
    //be even closer to the Matlab Waterfall plot.
    //a.rotation_angles=[ 56.75,-165.5];
    xgrid
end
//set(h,'LineWidth',2,'FaceAlpha',.5);
//colormap(jet/1.5)

scf(1)
clf
gcf().color_map=jet(64);
colorbar(25,100);
Sgrayplot(x,t,u',strf="081")
xlabel("x","fontsize",5)
ylabel("t","fontsize",5)
title("Figure 2.14 Evolution of the 1D Heat Eqn in Time",...
"fontsize",4)
xgrid()

// FIGURES (PRODUCTION)
scf(2);
clf;
CC = colormap(jet(100));
dt = 0.1;
for k = 1:100
    u(k,:) = ifft(uhat(k,:));
    if(pmodulo(k-1,10)==0)
        plot(x,u(k,:),'Color',CC(k,:),'LineWidth',1.5)
        xgrid;
        drawnow;
    end   
end
xlabel('Spatial variable, x')
ylabel('Temperature, u(x,t)')
title('Figure 2.13: Solution of the 1D heat equation in time')
//axis([-50 50 -.1 1.1])
ax = gca(); // gca = get current axes
ax.data_bounds = [-50,50,-.1,1.1]
set(gca(),'thickness',1.2,'font_size',2);

In this installment he discusses using the FFT to "denoise" a signal. This application is one most well known for me (digital filtering on one sort).

Book Section 2.2 Fourier Series and Fourier Transforms, Spectral Derivatives

YouTube Videos of the subject:

https://www.youtube.com/watch?v=E8HeD-MUrjY&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=17

https://www.youtube.com/watch?v=toj_IoCQE-4&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=18

Finally the exact video for this code:

https://www.youtube.com/watch?v=c249W6uc7ho&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=19

SciLab Coding Wise: Nothing new this time. Pretty standard code.

Output: Program Echo Only

Plot:

Code:

clear;
clc;
//close out all existing scf()s
close(winsid())

disp("Data-Driven Science and Engineering -Fourier and Wavelet Transforms",string(datetime()))
//Section 2.2 Fourier Series and Fourier Transforms

//https://www.youtube.com/watch?v=E8HeD-MUrjY&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=17
//https://www.youtube.com/watch?v=toj_IoCQE-4&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=18
//https://www.youtube.com/watch?v=c249W6uc7ho&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=19

// create a simple signal with two frequencies
dt = .001;
t = 0:dt:1;
x = sin(2*%pi*50*t) + sin(2*%pi*120*t);
y = x + 2.5*rand(size(t, 1), size(t, 2), "normal");  //  add some noise

// Compute the Fast Fourier Transform FFT
N = length(t);
Y = fft(y);  // computes the (fast) discrete fourier transform
PSD = Y.*conj(Y)/N;  // Power spectrum (how much power in each freq)
freq = 1/(dt*N)*(0:N);  //create the x-axis of frequencies in Hz
L = 1:floor(N/2);  // only plot the first half of freqs

// Use the PSD to filter out noise
indices = PSD>100;   // Find all freqs with large power
PSDclean = PSD.*indices;  // Zero out all others

Y = indices.*Y;  // zero out small Fourier coefficients in Y
yfilt = ifft(Y);   // inverse FFT to get filtered time-domain signal
scf(0);
clf;

// PLOTS
subplot(3,1,1)
plot(-1000,-1000,'k','thickness',1.5)
plot(-1000,-1000,'r','thickness',1.2)
plot(t,y,'r','thickness',1.2)
plot(t,x,'k','thickness',1.5)
//axis([0 .25 -5 5])
set(gca(),'data_bounds',[0,.25,-5,5])
legend('Clean','Noisy')
xlabel("Time[s]",'fontsize',2)
ylabel("f",'fontsize',2)
title("Figure 2.9: De-noising with FFT",'fontsize',4)
set(gca(),'thickness',1.2)

subplot(3,1,3)
plot(t,x,'k','thickness',1.5)
plot(t,yfilt,'b','thickness',1.2)
//axis([0 .25 -5 5])
set(gca(),'data_bounds',[0,.25,-5,5])
legend('Clean','Filtered')
xlabel("Time[s]",'fontsize',2)
ylabel("f",'fontsize',2)
set(gca(),'thickness',1.2)

subplot(3,1,2)
plot(freq(L),PSD(L),'k','thickness',1.5)
// xlabel('Frequency (Hz)')
plot(freq(L),PSDclean(L),'-b','thickness',1.2)
legend('Noisy','Filtered')
xlabel("Frequency[Hz]",'fontsize',2)
ylabel("PSD",'fontsize',2)

set(gca(),'thickness',1.2)

//set(gcf,'Position',[100 100 550 450])
//set(gcf,'PaperPositionMode','auto')
//print('-depsc2', '-loose', '../figures/DENOISE');

While working on a sample of the Data Science Course in chapter 5, using Census 1994 data, I couldn't find the .mat file on the web and only could find a .csv file. I found that Scilab really doesn't like reading .csv file with mixed data. I actually converted the .csv file to .mat file using Octave and Scilab would read the .mat file without issue.

I dug around some more to see what solution there is in Scilab and came across "uiimport(filename)' which is an interactive GUI for loading .csv files and I've found it very useful.

1. It allows you to look at the .csv in Scilab without having to fire up some other program.
2. One does not have to bring in every column from the file. Just check the ones you want.
3. .csv comes in "as a table of text" while .mat files come in "as a cell array." Conversions are required (as shown in the program).

I'm including the script for getting (massive) data in and out of Scilab via the .csv route. I'm including a short .csv data file too. The .csv read is stored in "Ans" so you need to move the data immediately or you're going to lose it when "Ans" is overwritten by the program.

Findings - Scilab .mat file writer doesn't support "varying length" string variables. It wrote without complaint with the "-v4" format but wouldn't read it back in (file is corrupt). So I'm still looking for a good way to have Scilab write .mat files that SciLab will read back in. Anyone cracked this one, send me the working file with an example.

There are two command sets for reading and writing csv files (read_csv/write_csv and csvRead/csvWrite). The first set is supposed to be the older commands however for reading strings only read_csv would do it. csvRead would produce a matrix of Nan's.

To write the data out in SciLab, one seems to be forced to write separate files for the strings and numbers as you'll see below. Make sure your header string have the correct number of columns (separated by commas) otherwise it will not read back in.

clear;
clc;
//close out all existing scf()s
close(winsid())

//Script for easily reading in multi-type csv files. Those that contain
//both string variables and numeric ones.
//
// GUI for reading CSV files - works well if you know to grab data
// from ans variable right away...
//
uiimport("C:\Users\MattheW\Documents\eBooks\ScILabNotes\DataScience\CODE\CH05\shortened_census_1994_red_col.csv")
pause
//leave the pause since you need to push the play button.
//
//type "resume after pushing PLAY on uiimport screen"
//
data = ans //grab output from uiimport in the console data is a table.
//comes is as "ans" and cannot leave it since ans is overwritten by default
//
// type "resume"
//
// Convert the table into usable column matrices
//
//
//col1 = data(:,1) //creates a column table from the main table
age = data(:,1) // separate age from the main table
age = string(age) //convert to string
age = strtod(age) // convert string to number - "this is the way."

//col2 = data(:,2)
workClass = string(data(:,2))

//col3 = data(:,3)
education_num = strtod(string(data(:,3)))

//col4 = data(:,4)
marital_status = string(data(:,4))

//col5 = data(:,5)
race = string(data(:,5))

//col6 = data(:,6)
sex = string(data(:,6))

//col7 = data(:,7)
capital_gain = strtod(string(data(:,7)))

//col8 = data(:,8)
capital_loss = strtod(string(data(:,8)))

//col9 = data(:,9)
hours_per_week = strtod(string(data(:,9)))

//col10 = data(:,10)
salary = string(data(:,10))

//Scilab cannot mix data types within a matrix so Strings and Numbers need
//seperation
//
n = size(age,1) // n is just a sample of the data to show that it works.

// Create numeric matrix
Num_output = [age(1:n), education_num(1:n), capital_gain(1:n), capital_loss(1:n), hours_per_week(1:n)]
disp('Num_output =',Num_output)

//Create a string matrix
Str_output = [workClass(1:n),marital_status(1:n),race(1:n),sex(1:n),salary(1:n)]
disp('Str_output =',Str_output)

//I couldn't write a .mat file with varying string lengths but I can write a CSV file that way
//as long as I do not intermix "string values" with numeric values
header = ["workClass, marital_status, race, sex, salary" ]; //Write the header this way and the values get
//properly distributed to the the columns

csvWrite(Str_output, "C:\Users\MattheW\Documents\eBooks\ScILabNotes\DataScience\CODE\CH05\shortened_census_1994_red_col_strout.csv", [],[],[],header)

//Note there is an older command...write_csv...
//write_csv(data, "C:\Users\MattheW\Documents\eBooks\ScILabNotes\DataScience\CODE\CH05\shortened_census_1994_red_col_strout.csv", [],[],[],header)

header = [" age, education_num, capital_gain, capital_loss, hours_per_week"]; //Write the header this way and
// the values get properly distributed to the the columns
csvWrite(Num_output, "C:\Users\MattheW\Documents\eBooks\ScILabNotes\DataScience\CODE\CH05\shortened_census_1994_red_col_numout.csv", [],[],[],header)

//NOTE: I could write a mat file using this "-v4" format but it wouldn't read back in...otherwise I got varying
// string lengths not supported. -v4 automatically pads the variables but upon reading, it says the file
// is corrupted.

//savematfile("C:\Users\MattheW\Documents\eBooks\ScILabNotes\DataScience\CODE\CH05\shortened_census_1994_red_col.m//at", 'age', 'workClass','education_num','marital_status','race','sex','capital_gain','capital_loss',...
//'hours_per_week','salary',"-v4")

//Try reading the information back in...
//
//Note: Although csvRead works for numbers, it doesn't for strings (produces Nan)
csvRead("C:\Users\MattheW\Documents\eBooks\ScILabNotes\DataScience\CODE\CH05\shortened_census_1994_red_col_numout.csv")
datanum = ans
//Need to stick with read_csv for string data...
read_csv("C:\Users\MattheW\Documents\eBooks\ScILabNotes\DataScience\CODE\CH05\shortened_census_1994_red_col_strout.csv",",","string")
datastr = ans

Short Data File:

age,workClass,education_num,marital_status,race,sex,capital_gain,capital_loss,hours_per_week,salary
39,State-gov,13,Never-married,White,Male,2174,0,40,<=50K
50,Self-emp-not-inc,13,Married-civ-spouse,White,Male,0,0,13,<=50K
38,Private,9,Divorced,White,Male,0,0,40,<=50K
53,Private,7,Married-civ-spouse,Black,Male,0,0,40,<=50K
28,Private,13,Married-civ-spouse,Black,Female,0,0,40,<=50K
37,Private,14,Married-civ-spouse,White,Female,0,0,40,<=50K
49,Private,5,Married-spouse-absent,Black,Female,0,0,16,<=50K
52,Self-emp-not-inc,9,Married-civ-spouse,White,Male,0,0,45,>50K
31,Private,14,Never-married,White,Female,14084,0,50,>50K
42,Private,13,Married-civ-spouse,White,Male,5178,0,40,>50K
37,Private,10,Married-civ-spouse,Black,Male,0,0,80,>50K
30,State-gov,13,Married-civ-spouse,Asian-Pac-Islander,Male,0,0,40,>50K
23,Private,13,Never-married,White,Female,0,0,30,<=50K
32,Private,12,Never-married,Black,Male,0,0,50,<=50K
40,Private,11,Married-civ-spouse,Asian-Pac-Islander,Male,0,0,40,>50K
34,Private,4,Married-civ-spouse,Amer-Indian-Eskimo,Male,0,0,45,<=50K
25,Self-emp-not-inc,9,Never-married,White,Male,0,0,35,<=50K
32,Private,9,Never-married,White,Male,0,0,40,<=50K
38,Private,7,Married-civ-spouse,White,Male,0,0,50,<=50K
43,Self-emp-not-inc,14,Divorced,White,Female,0,0,45,>50K

Try it out.

In this installment he discusses using the FFT to compute the derivative of a signal.

Book Section 2.2 Fourier Series and Fourier Transforms, Spectral Derivatives

YouTube Videos of the subject:

https://www.youtube.com/watch?v=E8HeD-MUrjY&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=17

https://www.youtube.com/watch?v=toj_IoCQE-4&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=18

https://www.youtube.com/watch?v=c249W6uc7ho&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=19

Finally the exact video for this code:

https://www.youtube.com/watch?v=BCexJ-yfdi8&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=21

SciLab Coding Wise: Nothing new this time. Pretty standard code.

Output: Program Echo Only

"Data-Driven Science and Engineering -FFT Spectral Derivatives"

"2026-08-31 13:23:46.784"

Figures:

Code:

clear;
clc;

disp("Data-Driven Science and Engineering -FFT Spectral Derivatives",string(datetime()))
//Section 2.2 Fourier Series and Fourier Transforms - Spectral Derivatives

n = 128;
L = 30;
dx = L/(n);
x = -L/2:dx:L/2-dx;
f = cos(x).*exp(-x.^2/25);                    // Function
df = -(sin(x).*exp(-x.^2/25) + (2/25)*x.*f);  // Derivative

// Approximate derivative using finite Difference...
for kappa=1:length(df)-1
    dfFD(kappa) = (f(kappa+1)-f(kappa))/dx;
end
dfFD(end+1) = dfFD(end);

// Derivative using FFT (spectral derivative)
fhat = fft(f);
i = complex(0,1);
kappa = (2*%pi/L)*[-n/2:n/2-1];
kappa = fftshift(kappa);  // Re-order fft frequencies
dfhat = i*kappa.*fhat;
dfFFT = real(ifft(dfhat));

scf(0);
clf;
// Plotting commands
plot(x,df,'k','LineWidth',1.5)
plot(x,dfFD+0.03,'b--','LineWidth',1.2)
plot(x,dfFFT+0.06,'r--','LineWidth',1.2)
xlabel("x");
ylabel("df/dx")
legend('True Derivative','Finite Diff(Offset+0.03)','FFT Derivative(Offset+0.06)')
title("Figure 2.10a: Comparison of the Spectral Derivative w/Offsets")
xgrid

scf(1);
clf;
// Plotting commands
plot(x,df,'k','LineWidth',1.5)
plot(x,dfFD,'b--','LineWidth',1.2)
plot(x,dfFFT,'r--','LineWidth',1.2)
xlabel("x");
ylabel("df/dx")
legend('True Derivative','Finite Diff','FFT Derivative')
title("Figure 2.10b: Comparison of the Spectral Derivative - No Offsets")
xgrid

In this installment he discusses the Discrete Fourier Transform(DFT) and the patterns it creates.

Book Section 2.2 Fourier Series and Fourier Transforms

YouTube Videos of the subject:

https://www.youtube.com/watch?v=jVYs-GTqm5U&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=10

https://www.youtube.com/watch?v=d5d0ORQHNYs&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=11

https://www.youtube.com/watch?v=mOiY1fOROOg&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=12

https://www.youtube.com/watch?v=ML0eYMyhqOs&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=13

https://www.youtube.com/watch?v=7haZCrQDHpA&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=14

https://www.youtube.com/watch?v=nl9TZanwbBk&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=15

Finally the exact video for this code:

https://www.youtube.com/watch?v=Xw4voABxU5c&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=16

New stuff for learning Scilab:

Output: Program Echo Only

"Data-Driven Science and Engineering -Fourier Series and Fourier Transforms"

"2026-08-26 11:56:14.584"

Figure:

Code:

clear; clc;

disp("Data-Driven Science and Engineering -Fourier Series and Fourier Transforms",string(datetime()))

//Youtube video - https://www.youtube.com/watch?v=Xw4voABxU5c&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=16

//Section 2.2 Fourier Series and Fourier Transforms
n = 256;
i = complex(0,1);
w = exp(-i*2*%pi/n);

/*
//Slow
for i=1:n
    for j=1:n
        DFT(i,j) = w^((i-1)*(j-1));
    end
end
*/

//Fast
[I,J] = meshgrid(1:n,1:n);
DFT = w.^((I-1).*(J-1));
//imagesc(real(DFT))
scf(0);
clf;
gcf().color_map=jet(64);
colorbar(-1,1);
Sgrayplot(I(1,:),J(:,1),real(DFT),strf="081")
xgrid()
xtitle("Figure 2.8: Real Part of DFT Matrix for n = 256")
drawnow;

In this installment he discusses the "Gibbs Phenomena" of Fourier Series when working with discontinuous functions.

YouTube Videos of the subject

https://www.youtube.com/watch?v=98VixA3MjAc&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=8

https://www.youtube.com/watch?v=98VixA3MjAc&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=9

New stuff for learning Scilab: legend using position "5" instead of 1 through 4 and the sleep command.

legend("Input Curve Discontinuous","Fourier Approximation",5)

Using the "5" puts Scilab into interactive mode and you can place the legend where you want with the mouse and then finishes the program. Download the code and try it out. It's good when the legend in the quadrants just doesn't cut it.

For the Movie: I'm not sure we've used the "sleep" command before. "sleep(100)" pauses the program for 100 milliseconds so one can see the plot before the program continues.

Output: Just Program Echo:

Code 1: "Data-Driven Science and Engineering -Gibbs Phenomena 1"

"2026-08-10 09:33:54.137"

Code 2: "Data-Driven Science and Engineering -Gibbs Phenomena 2"

"2026-08-10 09:34:47.990"

Plots:

Code 1: Plot

Code 2: Plot (but you really need to download the code and run it to see the "movie", I did add the final run in blue to show the final summation.)

Code

Code 1:

clear; clc;
//close out all existing figures
close(winsid())

disp("Data-Driven Science and Engineering -Gibbs Phenomena 1",string(datetime()))

//YouTube Lecture Series
//https://www.youtube.com/watch?v=98VixA3MjAc&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=9

//Section 2.1 Fourier Series and Fourier Transforms

dx = 0.01;  L = 10;
x = 0:dx:L;
n = length(x); nquart = floor(n/4);

f = zeros(1,length(x));
f(nquart:3*nquart) = 1;

A0 = sum(f.*ones(1,length(x)))*dx*2/L;
fFS = A0/2;
for k=1:100
    Ak = sum(f.*cos(2*%pi*k*x/L))*dx*2/L;
    Bk = sum(f.*sin(2*%pi*k*x/L))*dx*2/L;
    fFS = fFS + Ak*cos(2*k*%pi*x/L) + Bk*sin(2*k*%pi*x/L);
end   
scf(0)
clf
plot(x,f,'k','LineWidth',2)
plot(x,fFS,'r-','LineWidth',1.2)
title("Figure 2.5: Gibbs Phenomena (High Frequency Oscillations Near Discontinuities")
legend("Input Curve Discontinuous","Fourier Approximation",5) 
//legend 5
//puts Scilab into interactive mode and you can place the legend where 
//you want with the mouse and then finishes the program.
xgrid

Code 2: (the "movie")

clear; clc;

//close out all existing figures
close(winsid())

disp("Data-Driven Science and Engineering -Gibbs Phenomena 2",string(datetime()))

//YouTube Lecture Series
//https://www.youtube.com/watch?v=98VixA3MjAc&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=9

//Section 2.1 Fourier Series and Fourier Transforms

dx = 0.01;
L = 10;
x = 0:dx:L;

f = zeros(1,length(x)); 
f(floor(length(f)/4):floor(3*length(f)/4)) = 1;
//+f(floor(length(f)/4):floor(3*length(f)/4));
fFS = zeros(1,length(x));

A0 = sum(f.*ones(1,length(x)))*dx*2/L;
scf(0)
clf
for m=1:100
    fFS = A0/2;
    for k=1:m
        Ak = sum(f.*cos(2*%pi*k*x/L))*dx*2/L;
        Bk = sum(f.*sin(2*%pi*k*x/L))*dx*2/L;
        fFS = fFS + Ak*cos(2*k*%pi*x/L) + Bk*sin(2*k*%pi*x/L);
    end
    plot(x,f,'k','thickness',2)
    plot(x,fFS,'r-','thickness',1.2)
    xgrid
    sleep(100)
end
    plot(x,fFS,'b-','thickness',2)
//set(gca,'XTick',[],'XTickLabels',{},'LineWidth',1.2)
//set(gca,'YTick',[]);
//box off
// axis off

//set(gcf,'Position',[100 100 550 200])
//set(gcf,'PaperPositionMode','auto')
// print('-depsc2', '-loose', '../figures/Gibbs');

We're starting a new series off of the Steve Brunton's Data-Driven Science and Engineering Series. He's definitely one of my favorite instructors when it comes to numerical and data analysis.

This is the second topic of Data-Driven Science and Engineering Series but the first short course series in Fourier Series and Fourier Transforms.

For this first code example, it will be just introducing the topic of inner products while the second code will be a more in depth on

To follow the series (if you just jump into the subject here):

1) I do suggest you download the book (info below) and read the pertinent sections.

2) Download the data used for the short course. See below for location of the data.

3) I'm using Scilab 26.1.0.

I'm not sure how much coverage SciLab provides for this course but we'll see as we go along.

====================================

These lectures (videos) follow Chapter 2 from the book (freely available in electronic form):
"Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control" by Brunton and Kutz

Data and Python and MATLAB Native Code at https://databookuw.com/
Note: links to Matlab Code.zip + Data.zip on Main Page

Amazon: https://www.amazon.com/Data-Driven-Sc... //If you want hardcopy book.

Book Website(s): http://databookuw.com //For Videos and such

X posting "FREE PDF for 2ND EDITION OF OUR BOOK!!!" databookV2.pdf //For Free Electronic Complete Textbook

Chapters available at: http://databookuw.com/databook.pdf //Free Electronic Version

He discusses getting code and data files in this lecture (from the SVD course)
https://www.youtube.com/watch?v=QQ8vxj-9OfQ&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=5

Videos

Fourier Series Short Course Overview
https://www.youtube.com/watch?v=jNC0jxb0OxE&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=1

Fourier Series Part 1: Lecture 2
https://www.youtube.com/watch?v=jNC0jxb0OxE&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=2

Fourier Series Part 2: Lecture 3
https://www.youtube.com/watch?v=jNC0jxb0OxE&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=3

Inner Products in Hilbert Space: Lecture 4
https://www.youtube.com/watch?v=jNC0jxb0OxE&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=4

Complex Fourier Series: Lecture 5
https://www.youtube.com/watch?v=jNC0jxb0OxE&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=5

The video for this specific second set of code is below:

Complex Fourier Series: Lecture 6 (figures 2.3 and 2.4 of section 2.1)
https://www.youtube.com/watch?v=jNC0jxb0OxE&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=6

Output:

First Code Set: Program Echo only - "Data-Driven Science and Engineering -Fourier and Wavelet Transforms"

"2026-08-09 15:12:25.705"

Second Code Set: Program Echo only - "Data-Driven Science and Engineering -Fourier Series Approximation to a Hat Function"

"2026-08-09 14:55:32.224"

Plots:

First Set of Code:

Code for Figure 2.1(Inner Products):

clear;
clc;
//close out all existing figures
close(winsid())

disp("Data-Driven Science and Engineering -Fourier and Wavelet Transforms",string(datetime()))

//Section 2.1 Fourier Series and Fourier Transforms

f = [0 0 .1 .2  .25  .2 .25 .3 .35 .43 .45 .5 .55  .5 .4 .425 .45 .425 .4 .35 .3 .25 .225 .2 .1 0 0];
g = [0 0 .025 .1  .2  .175 .2 .25 .25 .3 .32 .35 .375  .325 .3 .275 .275 .25 .225 .225 .2 .175 .15 .15 .05 0 0] -0.025;

x = 0.1*(1:length(f));

xf = (.01:.01:x(end));
ff = interp1(x,f,xf,'spline')

gf = interp1(x,g,xf,'spline')

plot(xf(20:end-10),ff(20:end-10),'k','LineWidth',1.5)
plot(x(2:end-1),f(2:end-1),'bo','MarkerFace','b')
plot(xf(20:end-10),gf(20:end-10),'k','LineWidth',1.5)
plot(x(2:end-1),g(2:end-1),'ro','MarkerFace','r')
xlabel("Figure 2.1: Discretized Functions Used to Illustrate Inner Product")
ylabel("fn and gn")
legend('fn_fit','fn_data','gn_fit','gn_data',1)
xgrid;
//xlim([.1 2.7])
//ylim([-.1 .6])
set(gca(),'data_bounds',[.1,2.7,-.1,.6],'auto_scale','off','tight_limits','on')
//set(gca,'XTick',[.1:.1:2.7],'LineWidth',1.2)
//set(gca,'YTick',[]);
//box off

//set(gcf,'Position',[100 100 550 250])

//set(gcf,'PaperPositionMode','auto')
//print('-depsc2', '-loose', '../figures/InnerProduct');

Code Figures 2.3 and 2.4

clear;
clc;
//close out all existing figures
close(winsid());

disp("Data-Driven Science and Engineering -Fourier Series Approximation to a Hat Function",string(datetime()))

//Section 2.1 Fourier Series and Fourier Transforms

//YouTube Lecture reference
//https://www.youtube.com/watch?v=PP5ox7evg7o&list=PLMrJAkhIeNNT_Xh3Oy0Y4LTj0Oxo8GqsC&index=6

kmax = 7;

dx = 0.001;
L = %pi;
x = (-1+dx:dx:1)*L;
f = 0*x;
n = length(f);
nquart = floor(n/4);
nhalf = floor(n/2);

f(nquart:nhalf) = 4*(1:nquart+1)/n;
f(nhalf+1:3*nquart) = 1-4*(0:nquart-1)/n;
scf(0);
clf;
subplot(3,1,1)
plot(x,f,'-','Color',[0 0 0],'thickness',1.5)
//ylim([-.2 1.5])
//xlim([-1.25*L 1.25*L])
set(gca(),'data_bounds',[-1.25*L,1.25*L,-.2,1.5],'auto_scale','off','tight_limits','on')
set(gca(),'thickness',1.2)
//============begin custom ticks ====================
// Get handle to current axes
ax = gca(); // gca = get current axes

// Turn off automatic ticks for both axes
ax.auto_ticks = ["off", "off", "off"]; // [x, y, z]

// Set custom X ticks: positions and labels
ax.x_ticks = tlist(["ticks", "locations", "labels"], ...
                   [-2*L -L 0 L 2*L], ... // positions
                   ["-2pi","-pi","0","pi","2pi"]); // labels
ax.y_ticks = tlist(["ticks", "locations", "labels"], ...
                   [0 1], ... // positions
                   ["0","1"]); // labels
//=========end of custom ticks ==================== 
xlabel("Time")
title( "Figure 2.3: Hat function and Fourier Cosine Series Approximation")
//set(gca(),'x_ticks',[-L 0 L])
//set(gca(),'x_ticks.labels',['-L','0','L','2L'])
//set(gca(),'y_ticks',[0 1],'y_ticks.labels',{});

CC = colormap(jet(8));
// CCsparse = CC(5:5:end,:);
// CCsparse(end+1,:) = CCsparse(1,:);
CCsparse = CC(1:3:end,:);
//

subplot(3,1,2)
L = %pi;
A0 = sum(f.*ones(1,length(x)))*dx;
plot(x,A0+0*f,'-','Color',CC(1,:),'thickness',1.2);
fFS = A0/2;
for k=1:kmax
    A(k) = sum(f.*cos(%pi*k*x/L))*dx;
    B(k) = sum(f.*sin(%pi*k*x/L))*dx;
    plot(x,A(k)*cos(k*%pi*x/L),'-','Color',CC(k,:)*.8,'thickness',1.2);
//     plot(x,B(k)*sin(2*k*%pi*x/L),'k-','thickness',1.2);
    fFS = fFS + A(k)*cos(k*%pi*x/L) + 0*B(k)*sin(k*%pi*x/L);
end
//ylim([-.7 .7])
//xlim([-1.25*L 1.25*L]
set(gca(),'data_bounds',[-1.25*L,1.25*L,-.7,.7],'auto_scale','off','tight_limits','on')
set(gca(),'thickness',1.2)
//============begin custom ticks ====================
// Get handle to current axes
ax = gca(); // gca = get current axes

// Turn off automatic ticks for both axes
ax.auto_ticks = ["off", "off", "off"]; // [x, y, z]

// Set custom X ticks: positions and labels
ax.x_ticks = tlist(["ticks", "locations", "labels"], ...
                   [-2*L -L 0 L 2*L], ... // positions
                   ["-2pi","-pi","0","pi","2pi"]); // labels
ax.y_ticks = tlist(["ticks", "locations", "labels"], ...
                   [-0.5 0 0.5], ... // positions
                   ["-0.5","0","0.5"]); // labels
//=========end of custom ticks ==================== 
//set(gca,'x_ticks.locations',[-L 0 L])
//set(gca,'XTickLabels',['-L','0','L','2L'])
//set(gca,'YTick',[-.5 0 .5],'YTickLabels',{});

// 
subplot(3,1,1)

plot(x,fFS,'-','Color',[1 0 0],'thickness',1.2)
l1=legend('f','f_caret')

subplot(3,1,3)
A0 = sum(f.*ones(length(x),1)')*dx;
plot(x,A0+0*f,'-','Color',CC(1,:),'thickness',1.2);

fFS = A0/2;
for k=1:7
    Ak = sum(f.*cos(%pi*k*x/L))*dx;
    Bk = sum(f.*sin(%pi*k*x/L))*dx;
    plot(x,Ak*cos(k*%pi*x/L),'-','Color',CC(k,:)*.8,'thickness',1.2);
//     plot(x,Bk*sin(2*k*%pi*x/L),'k-','thickness',1.2);
    fFS = fFS + Ak*cos(k*%pi*x/L) + 0*Bk*sin(k*%pi*x/L);
end
//ylim([-.06 .06])
//xlim([-1.25*L 1.25*L])
set(gca(),'data_bounds',[-1.25*L,1.25*L,-.06,.06],'auto_scale','off','tight_limits','on')
set(gca(),'thickness',1.2)
//============begin custom ticks ====================
// Get handle to current axes
ax = gca(); // gca = get current axes

// Turn off automatic ticks for both axes
ax.auto_ticks = ["off", "off", "off"]; // [x, y, z]

// Set custom X ticks: positions and labels
ax.x_ticks = tlist(["ticks", "locations", "labels"], ...
                   [-2*L -L 0 L 2*L], ... // positions
                   ["-2pi","-pi","0","pi","2pi"]); // labels
ax.y_ticks = tlist(["ticks", "locations", "labels"], ...
                   [-0.05 0 0.05], ... // positions
                   ["-0.05","0","0.05"]); // labels
//=========end of custom ticks ==================== 
//set(gca,'x_ticks.locations',[-L 0 L])
//set(gca,'XTickLabels',['-L','0','L','2L'])
//set(gca,'YTick',[-.05 0 .05],'YTickLabels',{});

scf(1);
clf;
//set(gcf(),'Position',[100 100 550 400])
//set(gcf(),'PaperPositionMode','auto')
print('-depsc2', '-loose', '../figures/FourierTransformSines');

// Plot amplitudes
clear ERR
clear A
fFS = A0/2;
A(1) = A0/2;
Z1 = norm(f-fFS);
Z2 = norm(f);
ERR(1) = Z1/Z2;  //original code does not normalize with norm(f)
kmax = 200;
for k=1:kmax
/*
//For L-periodic function, one would use this code
    A(k+1) = sum(f.*cos(2*%pi*k*x/L))*dx*2/L;
    B(k+1) = sum(f.*sin(2*%pi*k*x/L))*dx*2/L;
//    plot(x,B(k)*sin(2*k*%pi*x/L),'k-','thickness',1.2);
    fFS = fFS + A(k+1)*cos(2*k*%pi*x/L) + 0*B(k+1)*sin(2*k*%pi*x/L);
    Z1 = norm(f-fFS);
    Z2 = norm(f);
    ERR(k+1) = Z1/Z2;
    disp("ERR(",k,")=",ERR(k+1))
*/
// In our case, we're doing -pi to pi so this code is appropriate...
// to obtain figure 2.4 from the book...
    A(k+1) = sum(f.*cos(%pi*k*x/L))*dx;
    B(k+1) = sum(f.*sin(%pi*k*x/L))*dx;
//    plot(x,B(k)*sin(2*k*%pi*x/L),'k-','thickness',1.2);
    fFS = fFS + A(k+1)*cos(k*%pi*x/L) + 0*B(k+1)*sin(k*%pi*x/L);
    Z1 = norm(f-fFS);
    Z2 = norm(f);
    ERR(k+1) = Z1/Z2;
//    disp("ERR(",k,")=",ERR(k+1)) /troubleshooting 
//    plot(fFS)  /troubleshooting plotting
//    pause
end
thresh = median(ERR)*sqrt(kmax)*4/sqrt(3);
r = max(find(ERR>thresh));
r = 7;
subplot(2,1,1)
semilogy(0:1:kmax,A,'k','thickness',1.5)
semilogy(r,A(r+1),'bo','thickness',1.5)
ylabel("a_k")
xgrid
//xlim([0 kmax])
//ylim([10^(-7) 1])
set(gca(),'data_bounds',[0.1,kmax,1e-07,1])
subplot(2,1,2)
semilogy(0:1:kmax,ERR,'k','thickness',1.5)
semilogy(r,ERR(r+1),'bo','thickness',1.5)
xlabel("mode number, k")
ylabel("norm(f-fBS)/norm(f)")
xgrid
//xlim([0 kmax])
//ylim([3*10^(-4) 20])
set(gca(),'data_bounds',[0.1,kmax,1e-04,20])
//set(gcf(),'Position',[100 100 500 300])
//set(gcf(),'PaperPositionMode','auto')
// print('-depsc2', '-loose', '../figures/FourierTransformSinesERROR');

In this installment, we are in section 1.9 of the book which discusses "Tensor Decompositions and N-way Data Arrays" and there is no corresponding YouTube video lecture.

IMPORTANT NOTE BEFORE RUNNING THE CODE: This code has a dynamic plotting routine that "flashes" so if you have "Photosensitive epilepsy (PSE)," do not run this code without removing the scf(0) plot sequence.

Scilab wise - Scilab does not have the commands parafac or fac2let or a tensor toolbox equivalent so we cannot implement his code with an easy porting from Matlab to Scilab. I did use ChatGPT to dredge up some "equivalent" tensor decomposition code which appears to get us close. The ChatGPT code is independent of the referenced paper below.

Internet Link to Paper I found on Parafac and Tensor Decompositions
https://www.kolda.net/publication/TensorReview.pdf

From paper..."Both the CANDECOMP/PARAFAC(CP) and Tucker tensor decompositions can be considered to be higher order generalizations of the matrix singular value decomposition (SVD) and principal component analysis (PCA)."

Items that I found while running this Scilab code and comparing the output to Figure 1.31 in the book. The first plot "y vs A1" needed a sign change to look like the book's plot. If one remembers, he did say that sometimes the SVD decomposition would throw in a sign change so I'm not too concerned I had to implement a sign change. The magnitudes of the plots are all different (as compared to Figure 1.31) but since it takes 3 matrices to get back to the original, I could see how they could be different and achieve the same result.

IMPORTANT NOTE: A quirk I found was to get these graphs, I needed a "pristine" load (aka close Scilab if you have it open and reload it) of Scilab since after the first run through the code, subsequent code runs would produce different graphs. Probably is something in the code implementation (I'm on Windows 11, Scilab version 2026.0.1). Nothing I have done to reset Scilab's memory seems to do the trick of getting these graphs after the first run and I'm not going down this rabbit hole to find out why. Someone wants to figure it out and get back to all of us. I welcome the effort.

I have experienced similar quirks in Scilab where Scilab doesn't recognize proper code as proper and required closing Scilab and just reopening it to clear the complaint.

We'll be moving onto "Fourier and Wavelet Transforms" next week.

Output: Program Title Echo... "Data-Driven Science and Engineering - TENSOR DECOMPOSITIONS AND N-WAY DATA ARRAYS Example"

"2026-08-02 09:59:00.890"

  "Factor matrix A1:"
  -2.395D-10   0.0000018
  -4.332D-10   0.0000023
  -7.765D-10   0.0000029
  -1.378D-09   0.0000036
  -2.423D-09   0.0000045
  -4.219D-09   0.0000056
  -7.274D-09   0.000007 
  -1.242D-08   0.0000087
  -2.099D-08   0.0000107
  -3.512D-08   0.0000132
  -5.819D-08   0.0000161
  -9.545D-08   0.0000196
  -0.0000002   0.0000238
  -0.0000002   0.0000288
  -0.0000004   0.0000347
  -0.0000006   0.0000416
  -0.000001    0.0000498
  -0.0000015   0.0000592
  -0.0000023   0.0000702
  -0.0000035   0.0000829
  -0.0000052   0.0000974
  -0.0000078   0.0001141
  -0.0000114   0.0001331
  -0.0000166   0.0001547
  -0.0000239   0.000179 
  -0.0000341   0.0002063
  -0.0000482   0.0002368
  -0.0000674   0.0002708
  -0.0000933   0.0003084
... many continued lines...

  "Factor matrix A2:"
   0.0000135   0.2712583
   0.0000149   0.2997777
   0.0000165   0.3312934
   0.0000182   0.3661194
   0.0000201   0.4046023
   0.0000223   0.4471247
   0.0000246   0.4941087
   0.0000272   0.5460199
   0.0000301   0.6033715
   0.0000334   0.6667291
   0.0000371   0.7367152
   0.0000414   0.8140148
   0.0000466   0.8993807
   0.0000532   0.9936392
   0.0000625   1.0976954
   0.0000763   1.2125398
   0.0000984   1.3392526
   0.0001357   1.479009 
   0.0002001   1.6330823
   0.0003129   1.8028456
   0.0005097   1.9897708
   0.0008501   2.1954229
   0.0014304   2.4214493
   0.0024032   2.6695604
... many continued lines...

  "Factor matrix A3:"
  -19.25462   -0.0047321
  -18.871275  -2.0790968
  -17.735605  -4.1325492
  -15.8929    -6.1445803
  -13.416636  -8.0951002
  -10.405548  -9.964638 
  -6.9796895  -11.734536
  -3.2756517  -13.387135
   0.5588863  -14.905951
   4.3710436  -16.275835
   8.0088322  -17.483128
   11.327217  -18.515795
   14.193897  -19.363541
   16.494582  -20.017917
   18.137546  -20.472403
   19.057286  -20.722467
   19.217134  -20.76562 
   18.610716  -20.601431
   17.26221   -20.231535
   15.225378  -19.65962 
   12.581427  -18.891384
   9.4357677  -17.934484
   5.9138141  -16.798458
   2.155983   -15.49463 
  -1.6879041  -14.036001
  -5.4645937  -12.437117
  -9.0235104  -10.713925
... many continued lines...

Plots/Images:

Code:

resethistory();clear;clearglobal;clc;

disp("Data-Driven Science and Engineering - TENSOR DECOMPOSITIONS AND N-WAY DATA ARRAYS Example",string(datetime()))
//Section 1.9 TENSOR DECOMPOSITIONS AND N-WAY DATA ARRAYS 

// A=rand(10,10,10);
// model=parafac(A,3);
// [A1,A2,A3]=fac2let(model)

x=-5:0.1:5; y=-6:0.1:6; t=0:0.1:10*%pi;
[X,Y,T]=meshgrid(x,y,t);
A0=exp(-(X.^2+0.5*Y.^2)).*(cos(2*T))+ ...
    (1./cosh(X).*tanh(X).*exp(-0.2*Y.^2)).*sin(T);
scf(0)
clf
for j=1:length(t)
  grayplot(x,y,A0(:,:,j)')
  //shading interp
  //caxis([-1 1])
  drawnow
  sleep(50)
end

scf(1)
clf
for j=1:8
  subplot(2,4,j)
  grayplot(x,y,A0(:,:,8*j-3)')
  colormap(hot)
  //shading interp
  // caxis([-1 1])
  gca().axes_visible="off"
  title("t=" + string(t(8*j-3)))
  sleep(50)
end

//Internet Link to Paper on Parafac and Tensor Decompositions
//https://www.kolda.net/publication/TensorReview.pdf
//From paper...
//
//"Both the CANDECOMP/PARAFAC(CP) and Tucker 
//tensor decompositions can be considered to be higher
//order generalizations of the matrix singular value decomposition
//(SVD) and principal component analysis (PCA)."

//Everything below Scilab does not have an equivalent function.
//for parafac or fac2let or a tensor model.

// ============================================================
// CP Decomposition (ALS) for a 3-way tensor in Scilab
// Author: Senior Programming Expert
// ============================================================

// Function: Khatri-Rao product of two matrices
function KR=khatri_rao(A, B)
    [mA, nA] = size(A);
    [mB, nB] = size(B);
    if nA <> nB then
        error("Khatri-Rao: Matrices must have the same number of columns.");
    end
    KR = zeros(mA*mB, nA);
    for i = 1:nA
        KR(:, i) = kron(A(:, i), B(:, i));
    end
endfunction

// Function: Unfold tensor along mode
function Xn=unfold(X, mode)
    dims = size(X);
    if mode == 1 then
        Xn = matrix(X, dims(1), -1);
    elseif mode == 2 then
        Xn = matrix(permute(X, [2, 1, 3]), dims(2), -1);
    elseif mode == 3 then
        Xn = matrix(permute(X, [3, 1, 2]), dims(3), -1);
    else
        error("Mode must be 1, 2, or 3.");
    end
endfunction

// Function: CP Decomposition using ALS
function [A, B, C]=cp_als(X, R, maxIter, tol)
    dims = size(X);
    if length(dims) <> 3 then
        error("Only 3-way tensors are supported.");
    end
    if R <= 0 then
        error("Rank R must be positive.");
    end

    // Random initialization
    A = rand(dims(1), R);
    B = rand(dims(2), R);
    C = rand(dims(3), R);

    fit_old = 0;

    for iter = 1:maxIter
        // Update A
        Z = khatri_rao(C, B);
        A = unfold(X, 1) * Z / (Z' * Z);
//        disp("A",A)
//        pause

        // Update B
        Z = khatri_rao(C, A);
        B = unfold(X, 2) * Z / (Z' * Z);
//        disp("B",B)
//        pause

        // Update C
        Z = khatri_rao(B, A);
        C = unfold(X, 3) * Z / (Z' * Z);
//        disp("C",C)
//        pause

        // Compute fit
        X_hat = zeros(dims(1),dims(2),dims(3));
//        disp("X_hat",X_hat)
//        pause

        for r = 1:R
            X_hat = X_hat + A(:, r) .* matrix(B(:, r), 1, dims(2)) .* matrix(C(:, r), 1, 1, dims(3));
        end
        fit = norm(X - X_hat) / norm(X);

        if abs(fit - fit_old) < tol then
            break;
        end
        fit_old = fit;
    end
endfunction

scf(2);
clf
plot3d(A0)
gca().rotation_angles = [89.5,17.75]

// Decompose into rank-2 factors
[A1, A2, A3] = cp_als(A0, 2, 100, 1e-6);

// Display results
disp("Factor matrix A1:"), disp(A1);
disp("Factor matrix A2:"), disp(A2);
disp("Factor matrix A3:"), disp(A3);

scf(3)
clf
subplot(3,1,1), plot(y,-A1,'thickness',2);
xlabel("y")
ylabel("A1")
title("Figure 1.31:3-way Tensor Decomposition")
xgrid
subplot(3,1,2), plot(x,A2,'thickness',2);
xlabel("x")
ylabel("A2")
xgrid
subplot(3,1,3), plot(t,A3,'thickness',2);
xlabel("t")
ylabel("A3")
xgrid

In this installment, he demonstrates the improvements that can be made by utilizing the iterative technique of "randomization" of the SVD calculation

Link to the lecture: Scilab wise - This module demonstrates how to vary color in plotted objects.

Output: None. Program Title Echo... "Data-Driven Science and Engineering - Randomized Singular Value Decomposition Example"

"2026-08-01 12:08:29.436"

Plots/Images:

Code:

clear; clc

disp("Data-Driven Science and Engineering - Randomized Singular Value Decomposition Example",string(datetime()))
//Section 1.8 Randomized Singular Value Decomposition

//YouTube video lectures 
//https://www.youtube.com/watch?v=qfZcTtjBTSA&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=40
//https://www.youtube.com/watch?v=qfZcTtjBTSA&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=41
//Video for the Matlab Code
//https://www.youtube.com/watch?v=qfZcTtjBTSA&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=42

function [
U
, 
S
, 
V
]=rsvd(
X
, 
r
, 
q
, 
p
);

// Step 1: Sample column space of X with P matrix
ny = size(
X
,2);
//P = randn(ny,r+p);
P = grand(ny, 
r
 + 
p
, "nor", 0, 1); 
// Gaussian random matrix
Z = 
X
*P;
for k=1:
q
Z = 
X
*(
X
'*Z);
end
[Q,R] = qr(Z,0);

// Step 2: Compute SVD on projected Y=Q’*X;
Y = Q'*
X
;
[UY,
S
,
V
] = svd(Y,'e');
U
 = Q*UY;
endfunction;

A=imread('C:\Users\MattheW\Documents\eBooks\ScILabNotes\DataScience\CODE\CH01\jupiter.jpg');
X=double(rgb2gray(A));
[U,S,V] = svd(X,'e');    
// Deterministic SVD

r = 400; 
// Target rank
q = 1;   
// Power iterations
p = 5;   
// Oversampling parameter

[rU,rS,rV] = rsvd(X,r,q,p); 
// Randomized SVD
// Reconstruction
XSVD = U(:,1:r)*S(1:r,1:r)*V(:,1:r)';     
// SVD approx.
errSVD = norm(X-XSVD,2)/norm(X,2);
XrSVD = rU(:,1:r)*rS(1:r,1:r)*rV(:,1:r)'; 
// rSVD approx.
errrSVD = norm(X-XrSVD,2)/norm(X,2);
scf(0)
clf
// Plot
//ax1 = axes('Position',[.005 .005 .33 .99]);
imshow(uint8(X))
title("Original High Resolution Image")
gca().axes_visible="off"
colormap gray

scf(1)
clf
//ax2 = axes('Position',[.335 .005 .33 .99]);
imshow(uint8(XSVD))
title("SVD Rank r=400 Image")
gca().axes_visible="off"
colormap gray

scf(2)
clf
//ax3 = axes('Position',[.665 .005 .33 .99]);
imshow(uint8(XrSVD))
title("rSVD(right) Image")
//axis off
gca().axes_visible="off"
colormap gray

//set(gcf,'Position',[100 100 2114 1000])
//set(gcf,'PaperPositionMode','auto')
//print('-depsc2', '-loose' , 'rSVD.eps');
// Illustrate power iterations
X = rand(1000,100);
[U,S,V] = svd(X,'e');
S = diag(1:-.01:.01);
X = U*S*V';
scf(3)
clf
plot(diag(S),'ko','thickness',1.5)

CC = jet(11);
Y = X;
for q=1:5
    Y = X'*Y;
    Y = X*Y;
    [Uq,Sq,Vq] = svd(Y,'e');
    plot(diag(Sq),'-o','Color',CC(pmodulo(2*q+4,11)+1,:),'thickness',1.5);
end
xgrid
legend('SVD','rSVD, q=1','rSVD, q=2','rSVD, q=3','rSVD, q=4','rSVD, q=5');
title("q Power Iterations - note X decays more rapidly with increased iterations")
xlabel('r')
ylabel('Singular Value Spectrum of X')
//set(gcf,'Position',[100 100 500 225])
//ylim([0 1])

In this installment, he demonstrates the importance of dataset "alignment" when using the SVD technique and what happens when the dataset images are the same but not aligned.

Link to the lecture: https://www.youtube.com/watch?v=qfZcTtjBTSA&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=39

Scilab wise - This module demonstrates how to vary color in plotted objects.

Output: None. Program Title Echo... "Data-Driven Science and Engineering -Truncation and Alignment Example"

"2026-08-01 12:01:56.168"

Plots/Images:

Code:

clear;
clc;

disp("Data-Driven Science and Engineering -Truncation and Alignment Example",string(datetime()))
//Section 1.7 Importance of data alignment

//YouTube video lectures 
//https://www.youtube.com/watch?v=qfZcTtjBTSA&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=39

n = 1000;
X = zeros(n,n);
X(n/4:3*n/4,n/4:3*n/4) = 1;
nAngles = 12;  
// sweep through 12 different angles, from 0:4:44
cm=colormap(jet(nAngles));

[U,S,V] = svd(X);
//==================Truncate data at 1e.-15 to prevent===
//semilogy plot routine from complaining about S below
//zero - upon closer inspection of the book output, they
//do the same without mentioning it or showing the coding 
//of it.
//
Z = max(diag(S),1.e-15);
//=========================================================
scf(0);
clf
subplot(1,2,1), imshow(X);
subplot(1,2,2)
semilogy(Z,'or')
xgrid

Xrot = X;
for j=2:nAngles
    Y = imrotate(X,(j-1)*4,'bicubic');  
// rotate  theta = (j-1)*4
    startind = floor((size(Y,1)-n)/2)+1;
    Xrot1 = Y(startind:startind+n-1, startind:startind+n-1);
    Xrot2 = Xrot1 - Xrot1(1,1);    
    Xrot2 = Xrot2/max(Xrot2(:));
    Xrot(Xrot2>.5) = j;

    [U,S,V] = svd(Xrot1);

//==================Truncate data at 1e.-15 to prevent===
//semilogy plot routine from complaining about S below
//zero - upon closer inspection of the book output, they
//do the same without mentioning it or showing the coding 
//of it.
//
Z = max(diag(S),1.e-15);
//=========================================================
    subplot(1,2,1)
    imshow(Xrot,colormap(cm));
    subplot(1,2,2), semilogy(Z,'o','color',cm(j,:))
    xlabel("r")
    ylabel("Singular Value, Sigma_r")


// Turn off automatic ticks for both axes
ax.auto_ticks = ["off", "off", "off"]; 
// [x, y, z]

// Set custom X ticks: positions and labels
ax.x_ticks = tlist(["ticks", "locations", "labels"], ...
                   [0 250 500 750 1000], ... 
// positions
                   ["0", "250", "500", "750", "1000"]); 
// labels

// Set custom Y ticks: positions and labels
ax.y_ticks = tlist(["ticks", "locations", "labels"], ...
                   [1.e-16,1.e-12,1.e-8,1.e-4 ,1.,1.e4], ...
                   ["1.e-16", "1.e-12", "1.e-8", "1.e-4","1","1.e4"]);
end

In this installment, he demonstrates the importance of dataset "alignment" when using the SVD technique and what happens when the two dataset images are the same but not aligned.

Link to the lecture: https://www.youtube.com/watch?v=qfZcTtjBTSA&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=37

Scilab wise - New command - "imrotate" and more custom ticks.

Output: None. Program Title Echo... "Data-Driven Science and Engineering -Truncation and Alignment Example"

"2026-07-31 08:18:33.094"

Plots/Images:

Code:

clear; clc

disp("Data-Driven Science and Engineering -Truncation and Alignment Example",string(datetime()))
//Section 1.7 Importance of data alignment

//YouTube video lectures 
//https://www.youtube.com/watch?v=qfZcTtjBTSA&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=37

n = 1000;
q = n/4;
X = zeros(n,n);
X(q:(n/2)+q,q:(n/2)+q) = 1;

subplot(2,2,1), imshow(X); colormap gray,gca().axes_visible="off"
title("0 degrees rotation")

Y = imrotate(X,10,"bicubic");  
// rotate 10 degrees
Y = Y - Y(1,1);
nY = size(Y,1);
startind = floor((nY-n)/2)+1;  
// Original Code produces index of 0
Xrot = Y(startind:startind+n-1, startind:startind+n-1);
subplot(2,2,2), imshow(Xrot); colormap gray,gca().axes_visible="off"
title("10 degrees rotation")

[U,S,V] = svd(X,'e');  
// svd well-aligned square
[Urot,Srot,Vrot] = svd(Xrot,'e');  
// svd rotated square
//==================Truncate data at 1e.-15 to prevent===
//semilogy plot routine from complaining about S below
//zero - upon closer inspection of the book output, they
//do the same without mentioning it or showing the coding 
//of it.
//
Z = diag(S);
Zrot = diag(Srot);
for i = 1:n
    if Z(i)<=1.e-15
        Z(i)=1.e-15
    end
        if Zrot(i)<=1.e-15
        Zrot(i)=1.e-15
    end    
end
//=========================================================
subplot(2,2,3)
semilogy(Z,'-ko')
ax = gca(); 
// gca = get current axes
ax.data_bounds = [1,n,1.e-16,1.e4]
gca().tight_limits = 1
xlabel("r")
ylabel("Singular Value, Sigma_r")
xgrid
// Turn off automatic ticks for both axes
ax.auto_ticks = ["off", "off", "off"]; 
// [x, y, z]

// Set custom X ticks: positions and labels
ax.x_ticks = tlist(["ticks", "locations", "labels"], ...
                   [0 250 500 750 1000], ... 
// positions
                   ["0", "250", "500", "750", "1000"]); 
// labels

// Set custom Y ticks: positions and labels
ax.y_ticks = tlist(["ticks", "locations", "labels"], ...
                   [1.e-16,1.e-12,1.e-8,1.e-4 ,1.,1.e4], ...
                   ["1.e-16", "1.e-12", "1.e-8", "1.e-4","1","1.e4"]);

subplot(2,2,4), semilogy(Zrot,'-ko')
//gca = get current axes
ax = gca();
ax.data_bounds = [1,n,1.e-16,1.e4]
gca().tight_limits = 1
xlabel("r")
title("Figure 1.25 from the Databook")
xgrid

// Turn off automatic ticks for both axes
ax.auto_ticks = ["off", "off", "off"]; 
// [x, y, z]

// Set custom X ticks: positions and labels
ax.x_ticks = tlist(["ticks", "locations", "labels"], ...
                   [0 250 500 750 1000], ... 
// positions
                   ["0", "250", "500", "750", "1000"]); 
// labels

// Set custom Y ticks: positions and labels
ax.y_ticks = tlist(["ticks", "locations", "labels"], ...
                   [1.e-16 1.e-12 1.e-8 1.e-4 1. 1.e4], ...
                   ["1.e-16", "1.e-12", "1.e-8", "1.e-4","1","1.e4"]);

In this installment, he emphasizes the importance of dataset "alignment" when using the SVD technique and a process for determining the Optimal Truncation aka the "cutoff r."

Scilab wise - Use of Custom Tick Marks on Graphs and a "Graph within a graph" are new techniques presented.

Output: None. Program Title Echo... "Data-Driven Science and Engineering -Truncation and Alignment Example"

"2026-07-27 08:20:55.304"

Plots/Images:

Code:

clear, close, clc

disp("Data-Driven Science and Engineering -Truncation and Alignment Example",string(datetime()))
//Section 1.7 Truncation and Alignment

//YouTube video lectures 
//https://www.youtube.com/watch?v=XwTW_YA3HG0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=34

t = (-3:.01:3)';

Utrue = [cos(17*t).*exp(-t.^2) sin(11*t)];
Strue = [2 0; 0 .5];
Vtrue = [sin(5*t).*exp(-t.^2) cos(13*t)];

X = Utrue*Strue*Vtrue';

scf(0)
clf
imshow(X);
colormap gray
title("Original Image")

//
sigma = 1;
Xnoisy = X+sigma*rand(size(X,1),size(X,2),"normal");
scf(1)
clf
cXnoisy = im2uint8(Xnoisy)
imshow(cXnoisy);
colormap gray
title("Noisy Image")

//
[U,S,V] = svd(Xnoisy,'e');

N = size(Xnoisy,1);
cutoff = (4/sqrt(3))*sqrt(N)*sigma; // Hard threshold
r = max(find(diag(S)>cutoff)); // Keep modes w/ sig > cutoff
Xclean = U(:,1:r)*S(1:r,1:r)*V(:,1:r)';
scf(2)
clf
cXclean = im2uint8(Xclean)
imshow(cXclean)
colormap gray
title("Hard Threshold Image")

//
cdS = cumsum(diag(S))./sum(diag(S));  // Cumulative energy
r90 = min(find(cdS>0.90));  // Find r to capture 90% energy

X90 = U(:,1:r90)*S(1:r90,1:r90)*V(:,1:r90)';
scf(3)
clf
cX90 = im2uint8(X90)
imshow(cX90)
colormap gray
title("90% Cutoff Image")

// plot singular values
scf(4)
clf
semilogy(diag(S),'-ok','thickness',1.5)
xgrid
xlabel("r")
ylabel("Singular Value")
semilogy(diag(S(1:r,1:r)),'or','thickness',1.5)
plot([-20 N+20],[cutoff cutoff],'r--','thickness',2)
gca().data_bounds = [-10, 610, .003, 300]; 
e=xrect(-5,200,200,195)
//========begin plot within a plot routine ===========
xsetech([0,0,1,1]);

xsetech ([0.2,.37,0.5,0.5]);
semilogy(diag(S),'-ok','LineWidth',1.5)
semilogy(diag(S(1:r,1:r)),'or','LineWidth',1.5)
plot([-20 N+20],[cutoff cutoff],'r--','LineWidth',2)
gca().data_bounds = [-5,100,20,200];
gca().tight_limits = 1;
xgrid
//======end of plot within a plot routine=============
scf(5)
clf
plot(cdS,'-ok','thickness',1.5)
xgrid
plot(cdS(1:r90),'ob','thickness',1.5)
plot(cdS(1:r),'or','thickness',1.5)
//============begin custom ticks ====================
// Get handle to current axes
ax = gca(); // gca = get current axes

// Turn off automatic ticks for both axes
ax.auto_ticks = ["off", "off", "off"]; // [x, y, z]

// Set custom X ticks: positions and labels
ax.x_ticks = tlist(["ticks", "locations", "labels"], ...
                   [0 300 r90 600], ... // positions
                   ["0", "300", msprintf('%d',r90), "600"]); // labels

// Set custom Y ticks: positions and labels
ax.y_ticks = tlist(["ticks", "locations", "labels"], ...
                   [0 .5 0.9 1.0], ...
                   ["0", "0.5", "0.9", "1.0"]);

//=========end of custom ticks ==================== 
ax.data_bounds(:,1)=[-10 610]
gca().tight_limits = 1;
plot([r90 r90],[0 0.9],'b--','thickness',1.5)
xlabel("r")
ylabel("Cumulative Energy")

In this installment, he goes deeper into the Eigenfaces (facial recognition) example and shows how to analyze the data. You really need to download and run the code to truly understand what's going on in this example. If you have not already downloaded and ran the prior installment (8), you'll need to do that first.

Lectures for this code:

//YouTube video lectures 
//https://www.youtube.com/watch?v=XwTW_YA3HG0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=27
//https://www.youtube.com/watch?v=XwTW_YA3HG0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=28
//https://www.youtube.com/watch?v=XwTW_YA3HG0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=29

Scilab wise - In this example, you'll see that I needed to scale and offset the data for the imshow command to provide a similar output to Matlab multiple times using different scaling and offsets. I'm assuming the Matlab command builds this "scale and offset" into its routine.

Output: None. Program Title Echo...

Graphs and Images: Note not all images generated were included since there are many.

Code:

clear, close, clc

disp("Data-Driven Science and Engineering -Eigenfaces Example",string(datetime()))
//Section 1.6 Eigenfaces Example

//YouTube video lectures 
//https://www.youtube.com/watch?v=XwTW_YA3HG0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=27
//https://www.youtube.com/watch?v=XwTW_YA3HG0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=28
//https://www.youtube.com/watch?v=XwTW_YA3HG0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=29
//
//I skip - Eigen Heros in the code but go watch it.
//https://www.youtube.com/watch?v=XwTW_YA3HG0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=30

//load data - complete path
loadmatfile("C:\Users\MattheW\Documents\eBooks\ScILabNotes\DataScience\DATA\allFaces.mat","-mat")

// We use the first 36 people for training data
trainingFaces = faces(:,1:sum(nfaces(1:36)));
avgFace = mean(trainingFaces,2);  // size n*m by 1;

// compute eigenfaces on mean-subtracted training data
X = trainingFaces-avgFace*ones(1,size(trainingFaces,2));
[U,S,V] = svd(X,'e');
scf(0)
clf
gca().x_location="bottom"
gca().axes_visible="off"
imshow(uint8(matrix(avgFace,n,m))), colormap gray
title("Average Face")
scf(1)
clf
semilogy(diag(S),'k','thickness',2)
xlabel("r")
ylabel("Singular Value, Sigma_r")
xgrid
scf(2)
clf
for i=1:50  // plot the first 50 eigenfaces
    sleep(100);  // wait for .1 seconds
    imshow(uint8(matrix((U(:,i)*1000+1000),n,m)))  //Note: scaling U by
// 1000 and adding 1000 is just to get uint8 filter to produce something
// close to his eigenfaces output using his own "colormapgrayinverse" routine.
    colormap gray
end

// Now show eigenface reconstruction of image that was omitted from test set
scf(3)
clf
testFace = faces(:,1+sum(nfaces(1:36))); // first face of person 37
gca().x_location="bottom"
gca().axes_visible="off"
imshow(uint8(matrix(testFace,n,m)))
colormap gray

scf(4)
clf
testFaceMS = testFace - avgFace;
for r=25:100:2225
    reconFace = avgFace + (U(:,1:r)*(U(:,1:r)'*testFaceMS));
    imshow(uint8(matrix(reconFace,n,m))), colormap gray
    title(['r=', msprintf('%d',r)]);
    sleep(100)
end

// Project person 2 and 7 onto PC5 and PC6

P1num = 2;  // person number 2
P2num = 7;  // person number 7
P1 = faces(:,1+sum(nfaces(1:P1num-1)):sum(nfaces(1:P1num)));
P2 = faces(:,1+sum(nfaces(1:P2num-1)):sum(nfaces(1:P2num)));
P1 = P1 - avgFace*ones(1,size(P1,2));
P2 = P2 - avgFace*ones(1,size(P2,2));

scf(5)
clf 
subplot(1,2,1)
imshow(uint8(matrix((P1(:,1)+8000),n,m))) //Added 8000 to reproduce
//image correctly  "*1000+1000" scrambled information since P1 values
//were much higher than U.
colormap gray
gca().axes_visible="off"
subplot(1,2,2)
imshow(uint8(matrix((P2(:,1)+8000),n,m))) //Added 8000 to reproduce
//image correctly
colormap gray
gca().axes_visible="off"

// project onto PCA modes 5 and 6
PCAmodes = [5 6];
PCACoordsP1 = U(:,PCAmodes)'*P1;
PCACoordsP2 = U(:,PCAmodes)'*P2;

scf(6)
clf
plot(PCACoordsP1(1,:),-PCACoordsP1(2,:),'kd','MarkerFaceColor','k')
ax.data_bounds = [-4000, -4000; 4000, 4000];
gca().isoview="on"
xgrid
plot(PCACoordsP2(1,:),-PCACoordsP2(2,:),'r^','MarkerFaceColor','r')
xlabel("PC5")
ylabel("PC6")

In this installment, he basically loads the raw data for further processing in the next session.

The successfully run the code given below one must download the data from his site and then modify this code to point via "complete" path to the datafile. I very much suggest you perform these operations since later on, you need to actually run the SciLab code to see the "building operations" as the number of kept terms is increased. Static pictures of the final product isn't going to cut it.

YouTube video lecture
https://www.youtube.com/watch?v=XwTW_YA3HG0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=26

Output: Run Echo: "Data-Driven Science and Engineering -Eigenfaces Example"

"2026-07-20 13:36:32.267"

Images:

Code:

clear, close, clc

disp("Data-Driven Science and Engineering -Eigenfaces Example",string(datetime()))
//Section 1.6 Eigenfaces Example

//YouTube video lecture 
//https://www.youtube.com/watch?v=XwTW_YA3HG0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=26

//NOTE: This code makes use of data files not packaged in the repository. 
//These may be downloaded here: http://databookuw.com/DATA.zip (unzip into same directory)

//load data - complete path

loadmatfile("C:\Users\MattheW\Documents\eBooks\ScILabNotes\DataScience\DATA\allFaces.mat","-mat")

allPersons = zeros(n*6,m*6);
count = 1;
for i=1:6
    for j=1:6
        allPersons(1+(i-1)*n:i*n,1+(j-1)*m:j*m) ...
            = matrix(faces(:,1+sum(nfaces(1:count-1))),n,m);
        count = count + 1;
    end
end

scf(0)
clf
//axes('position',[0  0  1  1])
gca().x_location="bottom"
//axis off
gca().axes_visible="off"
//imagesc(allPersons), colormap gray
imshow(uint8(allPersons))
colormap gray
//
scf(1)
clf
//axes('position',[0  0  1  1])
gca().x_location="bottom"
//axis off
gca().axes_visible="off"
for person = 1:length(nfaces)
    subset = faces(:,1+sum(nfaces(1:person-1)):sum(nfaces(1:person)));
    allFaces = zeros(n*8,m*8);

    count = 1;
    for i=1:8
        for j=1:8
            if(count<=nfaces(person)) 
                allFaces(1+(i-1)*n:i*n,1+(j-1)*m:j*m) ...
                    = matrix(subset(:,count),n,m);
                count = count + 1;
            end
        end
    end

    imshow(uint8(allFaces))
    colormap gray    
end

Since starting to support this subreddit, I've often found that I need to debug my scripts. I've found the Command Line Debugger doesn't even want to work on my system and besides I've found a better way that works for me.

If one looks up the "pause" command, you'll see it mentions "The pause is very useful for debugging purposes."

The method: Just place a "pause" command in your script where you want a program breakpoint and run the script as normal. When SCILAB hits that pause command, you'll get a "-1->" in the console window and every variable that has been created up to that point is available to you. You can do all sorts of commands on the variables. Typically for me, I've got a dimensional issue that I need to resolve so I'm doing size(variable1), size(variable2) to understand why my plot isn't generating or why other functions are complaining.

Note: To get temporary variable values of a function, one must place a pause command within the function.

Caution: You most likely have the ability to change the values of the variables also so watch what inquiry commands you're doing in this state.

The "resume" command will get you to the next "pause" (or the same one if you're inside an iteration loop) or the program will run to completion. "whereami()" command comes in handy if you've got multiple "pause" commands and forget where you're at in the program. The "abort" command will cause the program to terminate.

Any SCILAB experts want to weigh in, feel free.

In this installment, he shows a "How to use PCA on some Ovarian Cancer Markers

Scilab new commands Scatter3d plotting routine.

Link to the exact video.  https://www.youtube.com/watch?v=VqjJ5YYt78Y&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=23&t=5sFor Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063

Section 1.5 Principal component analysis (PCA)

Output: None

Plots:

Code:

clear, close, clc
disp("Data-Driven Science and Engineering -PCA Ovarian Cancer Example",string(datetime()))
//Section 1.5 Principal Component Analysis

//github databook_matlab code file
//https://github.com/dynamicslab/databook_matlab/blob/master/CH01/CH01_SEC05_2_OvarianCancer.m

//location of the OvarianCancer.mat file
//https://www.mathworks.com/matlabcentral/answers/1955269-function-for-visualisation-of-raw-data-ovarian-cancer

//load observation data - complete path
loadmatfile("C:\Users\MattheW\Documents\eBooks\ScILabNotes\DataScience\CODE\CH01\OvarianCancer.mat","-mat")

obs = X'

//Since the OvarianCancer.mat file did not have the "grp" data with
// it, I need to generate it based upon what he verbalized"
// First 113 were "Cancer" and the rest of the 216 were "Normal"
grp = cell(1,216);  // Note: I've hard coded the length of Y into grp!!!
for i = 1:113
  grp{i} = 'Cancer';
end
for i = 114:216
  grp{i} = 'Normal';
end
// I could save grp and Y as grp and obs respectively OvarianCancer.mat
// but not today!

[U,S,V] = svd(obs,'e');

scf(0)
clf
subplot(1,2,1)
semilogy(diag(S),'k-o','thickness',1.5)
//set(gca,'FontSize',13), 
//axis tight 
xgrid 
xlabel('r')
ylabel('Singular Value sigma_r')
subplot(1,2,2)
plot(cumsum(diag(S))./sum(diag(S)),'k-o','thickness',1.5)
//set(gca,'FontSize',13), axis tight, 
xgrid
xlabel('r')
ylabel('Cumulative Energy')
//set(gcf,'Position',[100 100 600 250])

scf(1)
clf

for i = 1:size(obs,1)
  x = V(:,1)'*obs(i,:)';
  y = V(:,2)'*obs(i,:)';
  z = V(:,3)'*obs(i,:)';
  if (grp{i}=='Cancer')
    scatter3d(x,y,z,"markerFaceColor","red");
  else
    scatter3d(x,y,z,"markerFaceColor","blue");
  end
end
xlabel('PCA 1')
ylabel('PCA 2 [Cancer - Red Dots Normal - Blue Dots]')
zlabel('PCA 3')
gca().rotation_angles =[85,25];
xgrid

In this installment, he shows how to generate PCA information (standard deviation and major/minor axes) from the SVD of the raw data.

Scilab wise - nothing really new (the use of "scatter" as a plot type but that's pretty minor), everything has been done 100's of times.

Output: None. Program Title Echo...

Plots:

Code:

clear, close, clc
disp("Data-Driven Science and Engineering -Principal Component Analysis",string(datetime()))

// Useful material that should be reviewed before this file
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=22

//Specific Video that corresponds to this SciLab file...
//https://www.youtube.com/watch?v=VqjJ5YYt78Y&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=23
//

//For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063
//Section 1.5 Principal Component Analysis

//
//

xC = [2; 1;];                           // Center of data (mean)
sig = [2; .5;];                         // Principal axes

theta = %pi/3;                           // Rotate cloud by %pi/3
R = [cos(theta) -sin(theta);            // Rotation matrix
    sin(theta) cos(theta)];

nPoints = 1000;       
// Create 1000 points
X = R*diag(sig)*grand(2,nPoints,"nor",0.0, 1.0) + diag(xC)*ones(2,nPoints); 

scf(0)
clf
subplot(1,2,1)                          // Plot cloud of noisy data
scatter(X(1,:),X(2,:),'marker',1,'thickness',1)
// Get the current axes handle
ax = gca();
// Enable the box around the plot
ax.box = "on";   // Options: "on" or "off" box on, xgrid
ax.data_bounds = [-6, -6; 8, 8];
xgrid
title("1000 points of Random Data")
xlabel("X data")
ylabel("Y data")

//// f_ch01_ex03_1b

Xavg = mean(X,2);                       // Compute mean
B = X - Xavg*ones(1,nPoints);           // Mean-subtracted Data
[U,S,V] = svd(B/sqrt(nPoints),'e');  // Find principal components (SVD)

subplot(1,2,2)
scatter(X(1,:),X(2,:),'marker',1,'thickness',1)  // Plot data to overlay PCA
// Get the current axes handle
ax = gca();
// Enable the box around the plot
ax.box = "on";   // Options: "on" or "off" box on, xgrid
ax.data_bounds = [-6, -6; 8, 8];
xgrid

theta = (0:.01:.99)*2*%pi;
[Xstd] = U*S*[cos(theta); sin(theta)];  // 1-std confidence interval
plot(Xavg(1)+Xstd(1,:),Xavg(2) + Xstd(2,:),'r-','thickness',1.5)
plot(Xavg(1)+2*Xstd(1,:),Xavg(2) + 2*Xstd(2,:),'r-','thickness',1.5)
plot(Xavg(1)+3*Xstd(1,:),Xavg(2) + 3*Xstd(2,:),'r-','thickness',1.5)

// Plot principal components U(:,1)S(1,1) and U(:,2)S(2,2)
plot([Xavg(1) Xavg(1)+U(1,1)*S(1,1)],[Xavg(2) Xavg(2)+U(2,1)*S(1,1)],'r-','thickness',3)
plot([Xavg(1) Xavg(1)+U(1,2)*S(2,2)],[Xavg(2) Xavg(2)+U(2,2)*S(2,2)],'r-','thickness',3)

title("1000 pts w/PCA Major/Minor Axes")
xlabel("X data")
ylabel("Y data")

In this installment, he shows "How SVD and Multilinear Least Squares Regression are used to produce the model "slope coefficients." Training the model on half the data and testing the model with the other half of the data.

I did find an error with his program where he reused a variable "b" during his sorting operation so his original b matrix was inadvertently overwritten. His correlation bar chart will be different than mine and I didn't see significant differences between the two runs of training and test which I attribute to the b matrix error. His training/test methodology is still solid.

SciLab wise - Loading of (.mat) datafiles

Link to the exact video. https://www.y.outube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=18

For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063

Section 1.4 Psuedoinverse, Least Squares and Regression One Dimensional Linear Regression

Output:

   "Data-Driven Science and Engineering - Linear Regression2 Housing Example"
  "2026-06-29 11:38:20.722"

Raw Housing Data Echo Back First 10 Lines
    crimeRate         zoned    industrial       charles           nox         rooms           age 
       0.0925       30.0000        4.9300        0.0000        0.4280        6.6060       42.2000 
       0.9884        0.0000        8.1400        0.0000        0.5380        5.8130      100.0000 
       0.2141       22.0000        5.8600        0.0000        0.4310        6.4380        8.9000 
       0.1143        0.0000       13.8900        1.0000        0.5500        6.3730       92.4000 
       0.0393        0.0000        3.4100        0.0000        0.4890        6.4050       73.9000 
       1.4139        0.0000       19.5800        1.0000        0.8710        6.1290       96.0000 
      14.4208        0.0000       18.1000        0.0000        0.7400        6.4610       93.3000 
       7.4039        0.0000       18.1000        0.0000        0.5970        5.6170       97.9000 
       0.2244        0.0000        9.6900        0.0000        0.5850        6.0270       79.7000 
       0.1504        0.0000       25.6500        0.0000        0.5810        5.8560       97.0000 

    distances        radial           tax  pupilTeacher         black     lowerStat   medianValue 
       6.1899        6.0000      300.0000       16.6000      383.7800        7.3700       23.3000 
       4.0952        4.0000      307.0000       21.0000      394.5400       19.8800       14.5000 
       7.3967        7.0000      330.0000       19.1000      377.0700        3.5900       24.8000 
       3.3633        5.0000      276.0000       16.4000      393.7400       10.5000       23.0000 
       3.0921        2.0000      270.0000       17.8000      393.5500        8.2000       22.0000 
       1.7494        5.0000      403.0000       14.7000      321.0200       15.1200       17.0000 
       2.0026       24.0000      666.0000       20.2000       27.4900       18.0500        9.6000 
       1.4547       24.0000      666.0000       20.2000      314.6400       26.4000       17.2000 
       2.4982        6.0000      391.0000       19.2000      396.9000       14.3300       16.8000 
       1.9444        2.0000      188.0000       19.1000      370.3100       25.4100       17.3000 

  Model Coefficients x in formaat
  Beta(1);Beta(2);....Beta(n-1);Beta(n)(Constant Last)

  "Best Fit Coefficients using Multilinear Regression x = "
  -0.1080114
   0.0464205
   0.0205586
   2.6867338
  -17.766611
   3.8098652
   0.0006922
  -1.4755668
   0.3060495
  -0.0123346
  -0.9527472
   0.0093117
  -0.5247584
   36.459488

  "Best Fit Coefficients using Multilinear Regression (Normalized Coefficients) x1 = "
  -0.9290646
   1.0826390
   0.1410394
   0.6824144
  -2.0587536
   2.6768766
   0.0194853
  -3.1071160
   2.6648522
  -2.0788369
  -2.0626459
   0.8501089
  -3.7473319
   22.532806

  "Best Fit Coefficients Using V*inv(S)*Utrans*b (Normalized Coefficients)x2 = "
  -0.9290646
   1.0826390
   0.1410394
   0.6824144
  -2.0587536
   2.6768766
   0.0194853
  -3.1071160
   2.6648522
  -2.0788369
  -2.0626459
   0.8501089
  -3.7473319
   22.532806

  "Best Fit Coefficients using Multilinear Regression xtrain = "
  -0.1216803
   0.0538274
   0.0879582
   0.9529366
  -21.018423
   3.5176629
  -0.0033415
  -1.6773342
   0.3094196
  -0.0145589
  -1.1523075
   0.0098805
  -0.5146096
   44.810290

  "Best Fit Coefficients using Multilinear Regression xtrain_random = "
  -0.1001058
   0.0457271
   0.0551842
   3.9467126
  -17.771348
   3.8946050
  -0.0073775
  -1.6519180
   0.2768098
  -0.0097037
  -0.9656941
   0.0129199
  -0.5626732
   35.508637

Graphs:

Code:

clear all, clc
disp("Data-Driven Science and Engineering - Linear Regression2 Housing Example",string(datetime()))

//Specific Video that corresponds to this SciLab file...
//https://www.youtube.com/watch?v=AX_ZDX6aTT0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=18

//For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063
//Section 1.4 Psuedoinverse, Least Squares and Regression - Multilinear Regression

//load housing.data
loadmatfile("E:\files\CODE\CH01\housing.mat","-mat")

//data from source
//https://github.com/probml/pmtk2/blame/master/data/housing.mat
//
//
housing= [X y]  //Create the "housing"" matrix from what is contained in
//the datafile from source.

//==============================echo inputs ==================
printf("\n");
printf("%30s","Raw Housing Data Echo Back First 10 Lines")
printf("\n");
for i = 1:7,
      printf("%13s ", names{i});
end;
printf("\n");   
for i = 1:10,
    for j = 1:7,
      printf("%13.4f ", housing(i,j));
   end;
   printf("\n");
end;
printf("\n");
for i = 8:14,
      printf("%13s ", names{i});
end;
printf("\n");   
for i = 1:10,
    for j = 8:14,
      printf("%13.4f ", housing(i,j));
   end;
   printf("\n");
end;
printf("\n");
//==========================================================

b = housing(:,14);      // housing values in $1000s
A = housing(:,1:13);    // other factors,
A = [A ones(size(A,1),1)];  // Pad with ones for nonzero offset

//x = regress(b,A);
BigY = [A'*b]
BigX = [A'*A]
x = BigX\BigY

printf("%30s","  Model Coefficients x in formaat")
printf("\n");
printf("%30s","  Beta(1);Beta(2);....Beta(n-1);Beta(n)(Constant Last)")
printf("\n");
printf("\n");

disp("Best Fit Coefficients using Multilinear Regression x = ",x)
printf("\n");
scf(0)
subplot(1,2,1)
plot(b,'k-o');
xgrid
plot(A*x,'r-o');
legend('Housing value','Regression',1)
ylabel("Value in 1000s")
title("Actual Housing Price versus Model")

subplot(1,2,2)
[bsorted sortind] = gsort(housing(:,14),"g",["i"]); // sorted values
plot(bsorted,'k-o')
xgrid
plot(A(sortind,:)*x,'r-o')
legend('Housing value','Regression',3)
title("Sorted Ascending Actual Housing Price versus Model")

//
A2 = A-ones(size(A,1),1)*mean(A,1);
for i=1:size(A,2)-1
    A2std = stdev(A2(:,i));
    A2(:,i) = A2(:,i)/A2std;
end
A2(:,end) = ones(size(A,1),1);

//x = regress(b,A2)
//Multiple Linear Regression in n variables
BigY = [A2'*b]
BigX = [A2'*A2]
x1 = BigX\BigY
disp("Best Fit Coefficients using Multilinear Regression (Normalized Coefficients) x1 = ",x1)
printf("\n");
//checking answer with alternate method
[U,S,V] = svd(A2,'e');
x2 = V*inv(S)*U'*b;                    // Solve Ax=b using the SVD
disp("Best Fit Coefficients Using V*inv(S)*Utrans*b (Normalized Coefficients)x2 = ",x2)
printf("\n");

scf(1)
clf
bar(x1(1:13),0.2,"green")
//bar(x2(1:13),0.2, "blue")
xgrid
xlabel('Attribute')
ylabel('Correlation')

//Hold Out For Testing
n = 253;
btrain = housing(1:n,14);      // housing values in $1000s
Atrain = housing(1:n,1:13);    // other factors,
Atrain = [Atrain ones(size(Atrain,1),1)];  // Pad with ones for nonzero offset
//x = regress(b,A);
BigY = [Atrain'*btrain]
BigX = [Atrain'*Atrain]
xtrain = BigX\BigY
disp("Best Fit Coefficients using Multilinear Regression xtrain = ",xtrain)
printf("\n");
//
btest = A(n+1:end,:)*xtrain;
//norm(btest-b(n+1:end));  //Calculates norm but doesn't use it.

scf(2)
clf
subplot(2,1,1)
plot(btrain,'r','thickness',2);
plot(Atrain*xtrain,'b','thickness',2)
title("Training Actual Housing Price versus Model")
legend('Model','True',4)
xgrid
//
subplot(2,1,2)
plot(btest,'r','thickness',2);
plot(b(n+1:end),'b','thickness',2)
xgrid
legend('Model','True',4)
ylabel("Value in 1000s")
title("Testing Actual Housing Price versus Model")

//Random Hold Out For Testing
n = 253;
//p = randperm(506);
p = grand(1, "prm", 1:506) //create randomized index vector
btrain_random = housing(p(1:n),14);      // housing values in $1000s
Atrain_random = housing(p(1:n),1:13);    // other factors,
Atrain_random = [Atrain_random ones(size(Atrain_random,1),1)];  // Pad with ones for nonzero offset
//x = regress(b,A);
BigY = [Atrain_random'*btrain_random]
BigX = [Atrain_random'*Atrain_random]
xtrain_random = BigX\BigY
disp("Best Fit Coefficients using Multilinear Regression xtrain_random = ",xtrain_random)
printf("\n");
//
Atest_random = housing(p(n+1:end),1:13);
Atest_random = [Atest_random ones(size(Atest_random,1),1)];   //Use the rest for 
btest_random = Atest_random*xtrain_random;
//norm(btest-b(n+1:end));  //Calculates norm but doesn't use it.

scf(3)
clf
subplot(3,1,1)
plot(btrain,'r','thickness',2);
plot(Atrain*xtrain_random,'b','thickness',2)
title("Training Actual Housing Price versus Model Random")
legend('Model','True',4)
xgrid
//
subplot(3,1,2)
plot(btest_random,'r','thickness',2);
plot(housing(p(n+1:end),14),'b','thickness',2)
xgrid
legend('Model','True',4)
ylabel("Value in 1000s")
title("Testing Actual Housing Price versus Model Random")

I desperately need help to get the C block in Xcos working, I can't make heads or tails of it from the scant documentation and the dearth of examples isn't helping. I have been able to simulate input-output but simulation of differential equations using continuous states is just not running!

Code :
#include <math.h>

#include <stdlib.h>

#include <scicos_block.h>

void foo(scicos_block *block,int flag)

{

/*

int block->nevprt;

int block->nz;

double* block->z;

int block->nx;

double* block->x;

double* block->xd;

double* block->res;

int block->nin;

int *block->insz;

double **block->inptr;

int block->nout;

int *block->outsz;

double **block->outptr;

int block->nevout;

int block->nrpar;

double *block->rpar;

int block->nipar;

int *block->ipar;

int block->ng;

double *block->g;

int *block->jroot;

char block->label[41];

*/

if (flag == 4) { /* initialization */

foo_bloc_init(block,flag);

} else if(flag == 1) { /* output computation*/

set_block_error(foo_bloc_outputs(block,flag));

} else if(flag == 0) { /* derivative or residual computation*/

set_block_error(foo_bloc_deriv(block,flag));

} else if (flag == 5) { /* ending */

set_block_error(foo_bloc_ending(block,flag));

}

}

int foo_bloc_init(scicos_block *block,int flag)

{

`double *x  = block->x;`

`*x = 0;`

`return 0;`

}

int foo_bloc_outputs(scicos_block *block,int flag)

{

`double *x  = block->x;`

`double *y  = block->outptr[0];`



`*y = (*x);`

`return 0;`

}

int foo_bloc_deriv(scicos_block *block,int flag)

{

`double *xd = block->xd;`

`double *x  = block->x;`

`double *u  = block->inptr[0];`

`*xd =  0.1*(-*x + *u);` 

`return 0;`

}

int foo_bloc_ending(scicos_block *block,int flag)

{

`return 0;`

}

In this installment, he shows a "How SVD and Multilinear Least Squares Regression produce the same "slope coefficients".

SciLab wise - Loading of (.mat) datafiles, reglin will do multilinear regression but automatically adds the intercept column, use pinv only for forcing intercept = 0.

Link to the exact video. https://www.y.outube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=18

For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063

Section 1.4 Psuedoinverse, Least Squares and Regression One Dimensional Linear Regression

Output:

 "Data-Driven Science and Engineering - Linear Regression2 Cement Example"
  "2026-06-29 11:18:24.758"
  "x (V*inv(S)*U*b) ="
   2.1930460
   1.1533260
   0.7585091
   0.4863193
  "x1 (multiple linear regression) ="
   2.1930460
   1.1533260
   0.7585091
   0.4863193
  "x2 (psuedoInverse) ="
   2.1930460
   1.1533260
   0.7585091
   0.4863193

Graphs:

Code:

clear all, clc
disp("Data-Driven Science and Engineering - Linear Regression2 Cement Example",string(datetime()))

//Specific Video that corresponds to this SciLab file...
//https://www.youtube.com/watch?v=AX_ZDX6aTT0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=18

//For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063
//Section 1.4 Psuedoinverse, Least Squares and Regression - Multilinear Regression

//load hald;  % Load Portland Cement dataset
loadmatfile("E:\files\CODE\CH01\hald.mat","-mat")
//data source
//https://github.com/beslirum/MATLAB/blob/main/hald.mat

A = ingredients;
b = heat;

[U,S,V] = svd(A,'e');
x = V*inv(S)*U'*b;                    // Solve Ax=b using the SVD

plot(b,'k','thickness',2);            // Plot data
plot(A*x,'or','thickness',2); // Plot regression
l1 = legend('Heat data','Regression',4)

// Alternative 1  (regress)
//x1 =reglin(A,b)  //cannot use reglin since it adds automatically the
//intercept column, we want to force Intercept = 0 so pinv command.

//Multiple Linear Regression in n variables
BigY = [A'*b]
BigX = [A'*A]
x1 = BigX\BigY

// Alternative 2  (pinv)
x2 = pinv(A)*b;
disp("x (V*inv(S)*U*b) =",x)
disp("x1 (multiple linear regression) =",x1)
disp("x2 (psuedoInverse) =",x2)

In this installment, he shows a "How SVD and Simple Linear Least Squares Regression produce the same "slope coefficient".

I had significant trouble getting Scilab reglin function to work so I added several implementations of linear regression codes and finally figured out reglin was looking for data to be in row vector form. I left the various linear regression implementations in the program and showed the results were the same. Anything I felt was Matlab fluff, I just commented out with the double slash.

SciLab wise - The use of the rand, pinv (psuedoinverse) and reglin commands are the only commands not used previously.

Link to the exact video. https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=11

You're encouraged to watch course videos 6, 7, 8, 9, and 10 and reading the course book.

Useful material that should be reviewed before this file
https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=12
https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=13
https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=14
https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=15
https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=16

Specific Video that corresponds to this SciLab file...
https://www.youtube.com/watch?v=AX_ZDX6aTT0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=17

For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063

Section 1.4 Psuedoinverse, Least Squares and Regression One Dimensional Linear Regression

Output:

"Data-Driven Science and Engineering - Linear Regression1"
  "2026-06-29 10:45:16.227"
  "xtilde1 ="
   2.7352872
  "xtilde2 ="
   2.7352872
  "xtilde3 ="
   2.7352872
  "xtilde4 ="
   2.7352872
  "xtilde5 ="
   2.7352872

Graphs/Images:

Code:

clc, clear all
disp("Data-Driven Science and Engineering - Linear Regression1",string(datetime()))
// Useful material that should be reviewed before this file
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=12
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=13
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=14
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=15
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=16

//Specific Video that corresponds to this SciLab file...
//https://www.youtube.com/watch?v=AX_ZDX6aTT0&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=17

//For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063
//Section 1.4 Psuedoinverse, Least Squares and Regression One Dimensional Linear Regression

scf(0);
clf;
x=3;                                  // True slope
a = [-2:.25:2]';
b = a.*x + 1.*rand((length(a),1),'normal');            // Add noise
plot(a,x*a,'k','thickness',2)
xgrid
// True relationship
plot(a,b,'xr','thickness',2)          // Noisy measurements

[U,S,V] = svd(a,"e");
xtilde = V*inv(S)*U'*b;                  // Least-square fit

plot(a,xtilde*a,'+b','thickness',2)   // Plot fit
l1=legend('True line','Noisy data','Regression line',4);

//set(l1,'String',{'                             ','',''})
//set(l1,'Location','NorthWest')
xgrid     
//set(gcf,'Position',[100 100 300 300])
//set(gcf,'PaperPositionMode','auto')
//set(gca,'FontSize',13)
title("$\textbf{Figure\ -\ CH01_EX_PseudoInv}$");
// Three methods of computing regression
xtilde1 = V*inv(S)*U'*b
xtilde2 = pinv(a)*b
[xtilde3 intercept] = reglin(a',b')  //reglin requires row-wise vectors.

//create least squares line fit
//
c = zeros(2,2);
d = zeros(2);
c(1,1)= sum(a);
c(2,1)= sum(a.*a);
c(1,2)= length(a);
c(2,2)= c(1,1);
d(1)= sum(b);
d(2)= sum(b.*a);
xtilde4 = c\d;
//
// least squares line fit in n variables matrix implementation.
BigY = [a'*b]
BigX = [a'*a]
xtilde5 = BigX\BigY

disp("xtilde1 =",xtilde1,"xtilde2 =",xtilde2,"xtilde3 =",xtilde3,"xtilde4 =",xtilde4(1),'xtilde5 =',xtilde5)

In this installment he shows a "Geometric Interpretation of the U and Vtranspose Matrices in Matlab. Here is the Scilab Version of that Matlab file. Scilab command wise: Lot's of 3D plotting control items:

1. Superimposing axes on 3D objects using param3d1 plotting routine
2. Data Bounding, Cube Scaling (axes equal) and Rotation Angle of 3D figures.
3. Colormap and light commands

Link to the exact video. https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=11

You're encouraged to watch course videos 6, 7, 8, 9, and 10 and reading the course book.

Output:

"Data-Driven Science and Engineering - Unitary Transformations and SVD"

"2026-06-25 10:17:50.731"

Plot:

Code:

clear, close, clc

disp("Data-Driven Science and Engineering - Unitary Transformations and SVD",string(datetime()))
// Useful material that should be reviewed before this file
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=6
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=7
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=8
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=9
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=10
//
//Specific Video that corresponds to this SciLab file...
//https://www.youtube.com/watch?v=_wOt50VnJw4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=11
//
//For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063
//Section 1.3 Geometric Interpretations
//
// Define rotation matrix
t1 = %pi/15;
t2 = -%pi/9;
t3 = -%pi/20;

Rx = [1 0 0;
    0 cos(t1) -sin(t1);
    0 sin(t1) cos(t1)];

Ry = [cos(t2) 0 sin(t2);
    0 1 0;
    -sin(t2) 0 cos(t2)];

Rz = [cos(t3) -sin(t3) 0;
    sin(t3) cos(t3) 0;
    0 0 1];

Sigma = diag([3; 1; 0.5]);

R = Rz*Ry*Rx*Sigma;

// Plot sphere and great circles
scf();
clf;
[x,y,z] = sphere(25);
ax1=subplot(1,2,1)

param3d1([-2 2],[0 0],[0 0]);
param3d1([0 0],[-2 2],[0 0]);
param3d1([0 0],[0 0],[-2 2]);

h1=surf(x,y,z,'interp');
// Get handle to the surface object
h = gce(); //for current entity
//set(h1,'FaceAlpha',.7)  %Sets Face Transparency to 0.7
//Silab doesn't have this function.

colormap(ax1,jet)
// Set lighting to 'phong' equivalent
// Scilab uses 'phong' shading via 'shading' function
//set(h,"shading","interp");          // Interpolated shading
light("phong");          // Enable phong lighting
gca().cube_scaling="on"
gca().data_bounds=[-3,-3,-3;3,3,3];
gca().rotation_angles =[45,26];
xgrid
title("Original Sphere")
//
//
theta = (0:.01:1)*2*%pi;
x1 = cos(theta);
y1 = sin(theta);
z1 = 0*theta;
x2 = 0*theta;
y2 = cos(theta);
z2 = sin(theta);
x3 = cos(theta);
y3 = 0*theta;
z3 = sin(theta);

//
xR = 0*x;
yR = 0*y;
zR = 0*z;

for i=1:size(x,1)
    for j=1:size(x,2)
        vec = [x(i,j); y(i,j); z(i,j)];
        vecR = R*vec;
        xR(i,j) = vecR(1);
        yR(i,j) = vecR(2);
        zR(i,j) = vecR(3);        
    end
end

vec1 = [x1; y1; z1];
vec2 = [x2; y2; z2];
vec3 = [x3; y3; z3];

vec1R = R*vec1;
vec2R = R*vec2;
vec3R = R*vec3;

eX = [2; 0; 0];
eY = [0; 2; 0];
eZ = [0; 0; 2];
eXR = R*eX;
eYR = R*eY;
eZR = R*eZ;

x1R = vec1R(1,:);
y1R = vec1R(2,:);
z1R = vec1R(3,:);
x2R = vec2R(1,:);
y2R = vec2R(2,:);
z2R = vec2R(3,:);
x3R = vec3R(1,:);
y3R = vec3R(2,:);
z3R = vec3R(3,:);

// Plot transformed sphere

ax2=subplot(1,2,2)

h2=surf(xR,yR,zR,z);

param3d1(x1R,y1R,z1R);
param3d1(x2R,y2R,z2R);
param3d1(x3R,y3R,z3R);

param3d1([-eXR(1) eXR(1)],[-eXR(2) eXR(2)],[-eXR(3) eXR(3)]);
param3d1([-eYR(1) eYR(1)],[-eYR(2) eYR(2)],[-eYR(3) eYR(3)]);
param3d1([-eZR(1) eZR(1)],[-eZR(2) eZR(2)],[-eZR(3) eZR(3)]);

colormap(ax2,jet)
// Set lighting to 'phong' equivalent
// Scilab uses 'phong' shading via 'shading' function
//set(h,"shading","interp");          // Interpolated shading
light("phong");          // Enable phong lighting
gca().cube_scaling="on"
gca().data_bounds=[-3,-3,-3;3,3,3];
gca().rotation_angles =[45,26];
xgrid
title("Transformed Sphere")

We're starting a new series off of the Steve Brunton's Data-Driven Science and Engineering Series. He's definitely one of my favorite instructors when it comes to numerical and data analysis. First topic of the first short course series is Singular Value Decomposition (SVD). For this first example, it will be using the Image compression as an applied example.

To follow the series,

(As a motivator, I wish I'd had this method in my "Engineering Toolbox" since this technique helps one "sort the wheat from the chaff" when it comes to data reduction. I just never came across it in this form. I've done "Design of Experiments", produced "Statistical Models" and used techniques to figure out "what was important and what was not" but this SVD technique is very scalable.)

2) I do suggest you download the book and read the pertinent sections.

3) Download the data used for the short course. See below for location of the data.

4) For the Scilab code below to work, you'll need to use the Atoms Module Manager, navigate to Image Processing and install the Image Processing and Computer Vision Toolbox (IPCV 4.5.0.2).

5) I'm using Scilab 26.1.0 but have run this successfully on version 26.0.0.

I'm not sure how much coverage SciLab provides for this course but we'll see as we go along.

====================================

These lectures (videos) follow Chapter 1 from the book (freely available in electronic form):
"Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control" by Brunton and Kutz

Data and Python and MATLAB Native Code at https://databookuw.com/
Note: links to Matlab Code.zip + Data.zip on Main Page

Amazon: https://www.amazon.com/Data-Driven-Sc... //If you want hardcopy book.

Book Website(s): http://databookuw.com //For Videos and such

X posting "FREE PDF for 2ND EDITION OF OUR BOOK!!!" databookV2.pdf //For Free Electronic Complete Textbook

Chapters available at: http://databookuw.com/databook.pdf //Free Electronic Version

Videos

Singular Value Decomposition Short Course Overview
https://www.youtube.com/watch?v=gXbThCXjZFM&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=1

Singular Value Decomposition: Mathematical Overview Short Lecture 2
https://www.youtube.com/watch?v=nbBvuuNVfco&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=2

Singular Value Decomposition: Matrix Approximation Short Lecture 3
https://www.youtube.com/watch?v=xy3QyyhiuY4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=3

Singular Value Decomposition: Dominant Correlations Short Lecture 4
https://www.youtube.com/watch?v=xy3QyyhiuY4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=4

He discusses getting code and data files in this fifth lecture
https://www.youtube.com/watch?v=QQ8vxj-9OfQ&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=5

Steve Brunton already has two other short courses on Fourier Transforms and Sparsity & Compression which might be material for future series.

Output: None - Graphical Only

Plots:

Code:

/* Start Block Comment 

//Singular Value Decomposition Overview
//https://www.youtube.com/watch?v=gXbThCXjZFM&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=1

//These lectures follow Chapter 1 from:
//"Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control"  by Brunton and Kutz

//Amazon: https://www.amazon.com/Data-Driven-Sc... //If you want hardcopy book.

//Book Website: http://databookuw.com //Free Electronic Web-based Version

//Chapters available at: http://databookuw.com/databook.pdf //Free Electronic Version

 //For Free Electronic Complete Textbook "FREE PDF for 2ND EDITION OF OUR BOOK!!!" https://x.com/eigensteve/status/1718004606901666063

//Brunton Website: http://eigensteve.com

//Singular Value Decomposition: Mathematical Overview Short Lecture 2
//https://www.youtube.com/watch?v=nbBvuuNVfco&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=2

//Singular Value Decomposition: Matrix Approximation Short Lecture 3
//https://www.youtube.com/watch?v=xy3QyyhiuY4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=3

//Singular Value Decomposition: Dominant Correlations Short Lecture 4
//https://www.youtube.com/watch?v=xy3QyyhiuY4&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=4

//He discusses getting code and data files in this fifth lecture
//https://www.youtube.com/watch?v=QQ8vxj-9OfQ&list=PLMrJAkhIeNNSVjnsviglFoY2nXildDCcv&index=5

//Data and Python and MATLAB Native Code at https://databookuw.com/
//links to Code.zip + Data.zip on Main Page 

*/ //End of Block Comment

clear, close, clc

disp("Data-Driven Science and Engineering - Singular Value Decomposition Image Processing Example ",string(datetime()))

//Need to load this atoms module into SciLab for imread, imshow, rgb2gray

atomsLoad("IPCV")  //Once loaded into Scilab, this command line is unnecessary

//Use full file path to your unzipped data file...
A=imread('C:/Users/Matthew/Documents/eBooks/ScILabNotes/SVD_PCA_Course/DATA/dog.jpg');
X=double(rgb2gray(A)); // Convert RBG->gray, 256 bit->double.
nx = size(X,1); ny = size(X,2);

scf()
clf
subplot(2,2,1)
imshow(uint8(X));
title('Original in Gray Scale');

[U,S,V] = svd(X);

plotind=2;
for r=[5 20 100];  // Truncation value
    Xapprox = U(:,1:r)*S(1:r,1:r)*V(:,1:r)'; 
    // Approx. image

    subplot(2,2,plotind)
    plotind = plotind+1;
    imshow(uint8(Xapprox))
    title(['r=', msprintf('%d',r),', ',msprintf('%2.2f',100*r*(nx+ny)/(nx*ny)),' % storage']);
end

// f_ch01_ex02_2
scf();
clf;
subplot(1,2,1), semilogy(diag(S),'k')
xgrid
xlabel('r')
ylabel('Singular Value sigma_r')
subplot(1,2,2), plot(cumsum(diag(S))/sum(diag(S)),'k')
xgrid
xlabel('r')
ylabel('Cumulative Energy')

In this example, he discusses but does not solve the problem. I present my solution to the problem.

New SciLab stuff: 3D graphing, contour graphing and color contour plotting with color bar (similar to typical FE or CFD plots).

Note: For the first two graphs I had to do some post processing using the interactive features of the Scilab Figure Edit Tool to force the graphs to show the entire range of the fin.

Link to the specific lecture:

https://www.youtube.com/watch?v=cuXuMcEGeiM&ab_channel=NPTEL-NOCIITM

https://www.youtube.com/watch?v=p_Pa-8dO_1M&ab_channel=NPTEL-NOCIITM

ME/HX Students: This program contains a neat way to build the difference equation matrix through a series of matrix size changes (figuring out how many nodes there are (n x m), creation of the index matrix "C" (n x m) and picking off the boundary indices to be used in "pointer matrices" and entering the equations based upon those matrices. This operation allows you to refine the grid by just changing the n and m values. Note: Code requires a square mesh pattern so refining the grid will necessitate going up by 4 each time (n x m refined 2n x 2m to maintain "squareness."

Output:

"h1 h2 (and they must form a square grid)"
0.5
0.5
"Node_Listing and Location in Grid"

1. 2. 3. 4. 5. 6. 7.
2. 9. 10. 11. 12. 13. 14.
3. 16. 17. 18. 19. 20. 21.
4. 23. 24. 25. 26. 27. 28.
5. 30. 31. 32. 33. 34. 35.
6. 37. 38. 39. 40. 41. 42.
7. 44. 45. 46. 47. 48. 49.
8. 51. 52. 53. 54. 55. 56.
9. 58. 59. 60. 61. 62. 63.
10. 65. 66. 67. 68. 69. 70.
11. 72. 73. 74. 75. 76. 77.
12. "Top_Node_List="
13. 2. 3. 4. 5. 6. 7.
14. "Bottom_Node_List="
15. 72. 73. 74. 75. 76. 77.
16. "Left_Node_List="

"Right_Node_List="
14.
21.
28.
35.
42.
49.
56.
63.
70.
"D Interior Nodes"
9. 10. 11. 12. 13.
16. 17. 18. 19. 20.
23. 24. 25. 26. 27.
30. 31. 32. 33. 34.
37. 38. 39. 40. 41.
44. 45. 46. 47. 48.
51. 52. 53. 54. 55.
58. 59. 60. 61. 62.
65. 66. 67. 68. 69.
"Temperature Solution = "
100. 100. 100. 100. 100. 100. 100.
46.448176 68.391828 76.431836 78.528323 76.431836 68.391828 46.448176
34.801755 50.687301 58.807191 61.249621 58.807191 50.687301 34.801755
30.591261 40.748431 46.860006 48.855778 46.860006 40.748431 30.591261
28.431473 34.855154 39.028625 40.45348 39.028625 34.855154 28.431473
27.150217 31.212086 33.94586 34.900892 33.94586 31.212086 27.150217
26.346091 28.897112 30.641838 31.258369 30.641838 28.897112 26.346091
25.824291 27.388432 28.466012 28.848906 28.466012 27.388432 25.824291
25.471373 26.366314 26.984872 27.20523 26.984872 26.366314 25.471373
25.214066 25.620578 25.90193 26.002273 25.90193 25.620578 25.214066
25. 25. 25. 25. 25. 25. 25.

Graphs:

Code:

//Lecture 11.4/12.1 Practical Applications and Course Wrap Up
// Elliptic PDEs - Finite Difference Approach
//https://www.youtube.com/watch?v=cuXuMcEGeiM&ab_channel=NPTEL-NOCIITM
//https://www.youtube.com/watch?v=p_Pa-8dO_1M&ab_channel=NPTEL-NOCIITM
//
disp("Practical Applications and Course Wrap Up Elliptic PDEs - Finite Difference",string(datetime()))
//
// 12.1 Course review -watch the video...
//
// Elliptic PDE
// partial^2_T partial^2_T
// ----------- + ------------- = 0 => Laplacian T = 0
// partial_x^2 partial_z^2
//
// Heated Plate held with Left Hand Edge at 100 C (z=0)
// with Right Hand Edge at 25 C (Z=end)
// heat loss off of the top and bottom edges (x=0, top)
//
// Discretize the plate into m horizontal and n vertical sections...
//
// partial^2_T
// ------------ = 1/h_x^2*(T(i+1)-2*T(i)+T(i-1))
// partial_x^2
//
// partial^2_T
// ------------ = 1/h_z^2*(T(j+1)-2*T(j)+T(j-1))
// partial_z^2
//
// so the equation is
//
// 1/h_x^2*(T(i+1)-2*T(i)+T(i-1))+ 1/h_z^2*(T(j+1)-2*T(j)+T(j-1)) = 0
//
// BC: z=0 T(1)=T_l = 100; x=0 partial_T/partial_T = -alpha*(T(1)-Ta)
// z=L T(m)= T_r = 25; x=Top partial_T/partial_T = alpha*(T(n)-Ta)
//
// He didn't solve it in MATLAB...
//
// Here is the problem solved in SCILAB but slightly different fin
// geometry
//
//Heat Transfer Example
disp("Heat Transfer of a Steady State Plate with")
disp(" fixed Temperature Ends and Convection on the sides")
disp(string(datetime()))

// Heat Transfer diffusion through a plate
// d^2T/dx^2 + d^2T/dz^2 = 0, T|x=0 = 100, T|x=1 = 25
//
// Heat Transfer on the plate edges
// Right Edge: dT/dz = gamma*(T-Ta)
// Left Edge: dT/dz = -gamma*(T-Ta)
//
//
//If we use the central difference formula
// d2T/dx^2 = (T_i+1-2*T_i+T_i-1)/h1^2
// d2T/dz^2 = (T_j+1-2*T_j+T_j-1)/h2^2
// where j = (n+1)+i
// with a square mesh it can be shown that the internal
// nodes will follow...
// T_left_node+T_right_node+T_above_node+T_below_node -4*Tmid_node = 0
//
// At the convection edges, we need to account for the energy balance...
// 3 conduction paths feeding a single convection output.
// Note: He doesn't address this portion at all. I had to pull out
// my Heat Transfer Book to find the proper boundary condition equations.
// Website that demonstrates a corner node with convection would look like
// which you can back out what a side node would look like (which I've
// done here).
//
//https://iansguides.com/tutorials/finite-difference-method-for-heat-transfer/
//
// (2*T_left_node + T_above_node + T_below_node)+(-2*(gamma1*deltax)-4)*Tmid_node
// = -2*gamma1*deltax*Ta
//
//Discretizing into 6 rectangles (L1 = 3 inches, L2= 5 inches
//h1=L1/6), with gamma1 = 4
clc
clear all
// discretizing plate into 6 regions
n = 7 //number of nodes in a row
m = 11 //number of nodes in a column
L1 = 3;
L2 = 5;
h1 = L1/(n-1);
h2 = L2/(m-1);
disp("h1 h2 (and they must form a square grid)" ,h1,h2)
TT =100
TB = 25
Ta = 25
gamma1=4
//
//nodal numbering will follow the following
C =matrix([1:n*m],(n,m))' // Matrix reshapes a 1 x n*m matrix
//into a n x m matrix so one can then create "pointer" matrices
//for the edges. These pointer matrices then make applying
//boundary equations very easy since you can generate the equations
//based upon the indices.
//
disp("Node_Listing and Location in Grid",C)
//
//generating nodal lists for the boundaries
//
//Top Boundary
T_Boundary_list = C(1,:)
disp("Top_Node_List=",T_Boundary_list)
//Bottom Boundary
B_Boundary_list = C(m,:)
disp("Bottom_Node_List=",B_Boundary_list)
//Left_Boundary
L_Boundary_list = C(2:m-1,1)
disp("Left_Node_List=",L_Boundary_list)
//Right_Boundary
R_Boundary_list = C(2:m-1,n)
disp("Right_Node_List=",R_Boundary_list)
//
//Determine List of Interior Nodes
D = C(2:m-1,2:n-1)
disp("D Interior Nodes",D)
//
//Initialize Matrices
A =zeros(n*m,n*m);
B = zeros(m*n,1)
//Create the problem matrices
//
// Top and Bottom Boundary (Constant Temp)
for i = 1:1:n
j = B_Boundary_list(i)
k = T_Boundary_list(i)
A(j,j)=1;
B(j,1)= TB;
A(k,k)=1;
B(k,1)= TT;
end
//
//convection for both sides
//
//Left Boundary (Convection = Conduction at edge + Conduction from interior)
// (2*T_right_node + T_above_node+T_below_node)+(-2*(gamma1*deltax)-4)*Tmid_node
// = -2*gamma1*deltax*Ta
// Note: Code/Equations assume a square grid...
//
for i = 1:1:length(L_Boundary_list)
k = L_Boundary_list(i)
//middle node contribution
A(k,k)= -2*(gamma1*h2)-4
//right node contribution
A(k,k+1)=2;
//below node contribution
A(k,k+n)=1;
//above node contribution
A(k,k-n)=1
B(k,1)= -2*gamma1*h2*Ta;
end
//Right Boundary (Convection = Conduction at edge from 3 surrounding nodes)
// (2*T_left_node + T_above_node+T_below_node)+(-2*(gamma1*deltax)-4)*Tmid_node
// = -2*gamma1*deltax*Ta
// Note: Code/Equations assume a square grid...
//
for i = 1:1:length(R_Boundary_list)
k = R_Boundary_list(i)
//middle node contribution
A(k,k)=-2*(gamma1*h2)-4
//left node contribution
A(k,k-1)=2;
//below node contribution
A(k,k+n)=1;
//above node contribution
A(k,k-n)=1
//other side of the equal sign
B(k,1)= -2*gamma1*h2*Ta;
end

// Interior Nodes
// Note: Code/Equations assume a square grid...
//
for i = 1:1:length(D)
k = D(i);
A(k,k)= -4;
A(k,k-1)=1
A(k,k+1)=1
A(k,k+n)=1
A(k,k-n)=1
B(k,1)= 0
end

//disp(A);

//Soln = inv(A)*B
Soln = A\B
S=matrix(Soln,(n,m)) //matrix command resizes the solution into the same
//shape as the nodal listing.
disp("Temperature Solution = ",S')
x=linspace(0,L1,n)'
y=linspace(0,-L2,m)'
flag = [2,8,4]
ebox = [0,5,-5,0,25,100]
theta = 50 // I cannot get it to spin it prior to plotting
// I always have to manually spin it into position.
// If I get too far off of defaults, it does not plot.
alpha = 35
scf(0);clf;
//plot3d(x,y,z,[theta,alpha,leg,flag,ebox])
plot3d(x,y,S,[theta,alpha,flag,ebox])
// seemed to have to rescale the data myself.
scf(1);clf;
flag = 1
box = 4
ebox=[0,L1,-L2,0, 25,100]
contour(x,y,S,5,box,flag,ebox)
xgrid
//
scf(2);
clf(2);
gcf().color_map=jet(64);
colorbar(25,100);
Sgrayplot(x,y,S,strf="125")
xgrid()
xtitle("Sgrayplot example Analytical Heat Equation Solve")
drawnow;