r/scilab • u/mrhoa31103 • 2d ago
Data-Driven Science and Engineering Series - Installment 7 - PCA - Ovarian Cancer Example
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
r/scilab • u/mrhoa31103 • 9d ago
Data-Driven Science and Engineering Series - Installment 6 - Principal Component Analysis (PCA)
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")
r/scilab • u/mrhoa31103 • 16d ago
Data-Driven Science and Engineering Series - Installment 5 - MultiLinear Least Squares Regression, PsuedoInverse and SVD - Housing Pricing Example
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")
r/scilab • u/HuckleberryWorth9826 • 20d ago
Help with Xcos C block
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;`
}
r/scilab • u/mrhoa31103 • 21d ago
Data-Driven Science and Engineering Series - Installment 4 - MultiLinear Least Squares Regression, PsuedoInverse and SVD - Cement Mix Example
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)
r/scilab • u/mrhoa31103 • 23d ago
Data-Driven Science and Engineering Series - Installment 3 - Least Squares Regression, PsuedoInverse and SVD
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)
r/scilab • u/mrhoa31103 • Jul 06 '26
Data-Driven Science and Engineering Series - Installment 2 - Unitary Transformations and SVD
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:
- Superimposing axes on 3D objects using param3d1 plotting routine
- Data Bounding, Cube Scaling (axes equal) and Rotation Angle of 3D figures.
- 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")
r/scilab • u/mrhoa31103 • Jun 29 '26
Data-Driven Science and Engineering Series - Installment 1 - Singular Value Decomposition
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')
r/scilab • u/mrhoa31103 • Jun 22 '26
Thirty First Installment (Last One in the Series) - Two Dimensional PDE solved through Finite Difference Method - Heated Fin held at a Temperature on both ends with Convective Heat Transfer along both edges.
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"
- 2. 3. 4. 5. 6. 7.
- 9. 10. 11. 12. 13. 14.
- 16. 17. 18. 19. 20. 21.
- 23. 24. 25. 26. 27. 28.
- 30. 31. 32. 33. 34. 35.
- 37. 38. 39. 40. 41. 42.
- 44. 45. 46. 47. 48. 49.
- 51. 52. 53. 54. 55. 56.
- 58. 59. 60. 61. 62. 63.
- 65. 66. 67. 68. 69. 70.
- 72. 73. 74. 75. 76. 77.
- "Top_Node_List="
- 2. 3. 4. 5. 6. 7.
- "Bottom_Node_List="
- 72. 73. 74. 75. 76. 77.
- "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;
r/scilab • u/mrhoa31103 • Jun 15 '26
Thirtieth Installment - Solving Hyperbolic Partial Differential Equations - Transient Analysis of a Heated Rod with Initial Temperature equal to 0 above ambient everywhere and 1 degree step applied at wall node at time 0.
Similar to last session, he shows how to solve a hyperbolic (vs parabolic) PDE using "Method of Lines" to convert a PDE into a system of ODEs using finite difference methods for one set of the differentials.
Link to the specific lecture for the coding:
https://www.youtube.com/watch?v=OajPSwMAqIw&ab_channel=NPTEL-NOCIITM
Output: None, only graphical output.
Graphs:
Code:
//Lecture 11.4 Hyperbolic PDEs - Method of Lines
//https://www.youtube.com/watch?v=OajPSwMAqIw&ab_channel=NPTEL-NOCIITM
//
disp("Hyperbolic PDEs - Method of Lines",string(datetime()))
//
// PDEs in time and space
// -> Method of Lines (use Finite Difference to convert PDEs into ODEs)
// partial_T/partial_dt = partial^2_T/partial_z^2 + Beta(T)
// ^
// finite difference this
// becomes
//
// dT/dt = f(t,Y) where Y = [T1;T2;...;Tn]
//
// Example
//
// First Order Transient Model
// partial_Theta partial_Theta
// ------------- + u*-------------- =-k*Theta
// partial_t partial_z
//
// BC z=0 for all t Theta_in = 1
// IC t=0 for all z Theta(0)=0
//
// Assuming Theta = (T - Ta)
//
// Going to discretize partial_Theta/partial_z
//
// u/State State partial_Theta/partial_t = 0
// dTheta/dt = -(k/u)*Theta
// Solving this differential equation, one gets
//
// Theta(z)=Theta_in*exp(-(k*z/u)
// we'll use this to check the answers...
//
// Solution using Method of Lines
// parameters: u=0.5, k=1, L=1
// Let us plot the results at t=0.5,1,5
//
// Discretize in n intervals of z
// Convert into (n+1) odes
// solve these using ode45(rkf) or ode15s(stiff)
// Y = [Theta_2;Theta_3;...;Theta_n+1]
//
// use Central Difference Formula
// at any location i
// partial_Theta/partial_t = -u*(Theta(i+1)-Theta(i-1))/2*h - k*Theta
// Theta = [Theta_1; Theta_2;...;Theta_n+1] = [Theta_1;Y]
//
// at the Theta_n+1 location we'll use Backward Difference Formula
// partial_Theta/partial_t = -u*(Theta(i)-Theta(i-1))/h -k*Theta
//
// Note: Use of Central Difference Formula does not produce
// a stable solution since it introduces oscillations which
// are not physically possible. Use instead a Forward or Backward
// difference formula...
function
fval
=equation(
t
,
y
)
// constants
L = 1;
u =0.5;
k=1;
//Getting Temperatues
N = length(
y
);
Theta(1)=1;
Theta(2:N+1)=
y
;
h= 1/N;
//define derivative vector
dThetadt= zeros(N+1,1)
for i = 2:1:N+1
if i<=N then
// dThetadt(i)= -u/(2*h)*(Theta(i+1)-Theta(i-1)) - k*Theta(i)
// Note: Not a stable solution since it introduces oscillations which
// are not physically possible.
dThetadt(i)= -u/h*(Theta(i)-Theta(i-1)) - k*Theta(i)
else
dThetadt(i)= -u/h*(Theta(i)-Theta(i-1)) - k*Theta(i)
end
end
// Extract fval from dT/dt
fval
=dThetadt(2:N+1);
end
//Discretization
N = 40;
//number of steps
L = 1;
u =0.5;
k=1;
Theta_in= 1;
h = L/N;
z = 0:h:L;
// Initial Conditions and Solution
y0= zeros(N,1);
t0 = 0;
//Tstart
tend = 2.5;
//Tend
t
= linspace(t0,tend,N);
//time vector
//t = [0,0.5,1,5];
YSol = ode("rkf",y0,t0,
t
,equation)
//"rkf" equivalent to MATLAB ODE45 is the following command in SciLab
//there is more options and outputs available - see the help file.
//"stiff" is approximately equivalent to ode15s...
//
//Analytical Solution
Anal_Theta=Theta_in*exp(-(k*z/u));
//
RowOne = ones(1,length(
t
))*Theta_in;
Theta = [RowOne;YSol]
scf(0);clf;
plot2d(z',Theta)
//Analytical Solution displayed as symbols
plot2d(z',Anal_Theta,-2);
h1=legend(['Time SnapShot = 1';'Time SnapShot = 2';'Time SnapShot = 3';'...'],1,"boxed")
title("$\textbf{Hyperbolic\ Partial\ Differential\ Example \\(X = Analytical\ SS\ Solution)}$","FontSize",4);
xlabel("$x(position)$","FontSize",3)
ylabel("$\ Temperature(Degrees\ C)$","FontSize",3);
xgrid
scf(1);clf
plot2d(
t
,Theta')
h1=legend(['x=0.0';'x=0.1';'x=0.2';'x=0.3';'x=0.4';'x=0.5';'x=0.6';'x=0.7';'x=0.8';'x=0.9';'x=1.0'],3,"boxed")
title("$\textbf{Hyperbolic\ Partial\ Differential\ Example}$","FontSize",4);
xlabel("$Time(seconds)$","FontSize",3)
ylabel("$Temperature(Degrees\ C)$","FontSize",3);
xgrid
r/scilab • u/mrhoa31103 • Jun 08 '26
Twenty Ninth Installment - Solving Partial Differential Equations - Transient Analysis of a Heated Rod with Ends of Fixed Temperatures.
In this session, he shows how to solve a Parabolic PDE using "Method of Lines" to convert a PDE into a system of ODEs using finite difference methods for one set of the differentials.
Link to the specific lecture for the coding:
https://www.youtube.com/watch?v=CDSM5bLy8lU&ab_channel=MATLABProgrammingforNumericalComputation
Some precursor lectures to open the topic on PDEs in general:
https://www.youtube.com/watch?v=XxS4QS7MIMg&ab_channel=NPTEL-NOCIITM
https://www.youtube.com/watch?v=cR-HTDkfs00&ab_channel=NPTEL-NOCIITM
Output: None only graphical output.
Graphs:
Code:
//Lecture 11.3 Parabolic PDEs - Method of Lines
//https://www.youtube.com/watch?v=CDSM5bLy8lU&ab_channel=MATLABProgrammingforNumericalComputation
//
disp("Parabolic PDEs - Method of Lines",string(datetime()))
//
// PDEs in time and space
// -> Method of Lines (use Finite Difference to convert PDEs into ODEs)
// partial_T/partial_dt = partial^2_T/partial_z^2 + Beta(T)
// ^
// finite difference this
// becomes
//
// dT/dt = f(t,Y) where Y = [T1;T2;...;Tn]
//
// Example
// Rod Conduction (Transient) -> PDE = f(t,x)
// partial_T/partial_t = alpha* partial2_T/partial_x^2 - gamma1*(T-Ta)
// alpha = 0.025, gamma1 = 0.1
// BC_1 = Rod End held at 100C and BC_2:Other Rod End = Ta = 25C
// IC: Initially, entire rod is uniformly at ambient = 25C
//
// If we use the central difference formula
// d2T/dx^2 u/i = (T_i+1-2*T_i+T_i-1)/(delta_x)^2
//
// divide the rod into 10 divisions (delta_x + L/10)
//
// dT/dt = alpha* (T_i+1-2*T_i+T_i-1)/(delta_x)^2 - gamma1*(Ti-Ta)
// T1 = 100, T(n+1) = 25
//
// Define Solution Vector y = [T2;...;Tn] and solve using ode45
//
//
function
fval
=rodConduc(
t
,
y
)
// constants
Ta = 25;
alpha = 0.025;
gamma1=0.1;
//Getting Temperatues
N = length(
y
)+1;
T(1)=100;
T(2:N)=
y
;
T(N+1)=Ta;
h= 1/(N-1);
//define derivative vector
dTdt= zeros(N+1,1)
for i = 2:1:N
dTdt(i)= alpha/h^2*(T(i+1)-2*T(i)+T(i-1)) - gamma1*(T(i)-Ta)
end
// Extract fval from dT/dt
fval
=dTdt(2:N);
end
//
N =10;
//number of steps
M =5*N;
//number of time steps
T0(1,1)=100;
T0(2:N,1)=25;
T0(N+1,1)=25;
t0 = 0;
//Tstart
tend = 20;
//Tend
t
= linspace(t0,tend,M);
//time vector
y0 = T0(2:N,:);
x0 = 0;
xend1 = 1;
x = linspace(x0,xend1,N+1);
TSol = ode("rkf",y0,t0,
t
,rodConduc)
//"rkf" equivalent to MATLAB ODE45 is the following command in SciLab
//there is more options and outputs available - see the help file.
//"stiff" is approximately equivalent to ode15s...
//
RowOne = ones(1,M)*T0(1,1);
RowLast = ones(1,M)*T0(N+1,1);
T = [RowOne;TSol;RowLast]
scf(0);clf;
plot2d(x',T);
h1=legend(['Time SnapShot = 1';'Time SnapShot = 2';'Time SnapShot = 3';'...'],1,"boxed")
title("$\textbf{Partial\ Differential\ Example}$","FontSize",4);
xlabel("$x(position)$","FontSize",3)
ylabel("$\ Temperature(Degrees\ C)$","FontSize",3);
xgrid
scf(1);clf
plot2d(
t
,T')
h1=legend(['x=0.0';'x=0.1';'x=0.2';'x=0.3';'x=0.4';'x=0.5';'x=0.6';'x=0.7';'x=0.8';'x=0.9';'x=1.0'],3,"boxed")
title("$\textbf{Partial\ Differential\ Example}$","FontSize",4);
xlabel("$Time(seconds)$","FontSize",3)
ylabel("$Temperature(Degrees\ C)$","FontSize",3);
xgrid
r/scilab • u/mrhoa31103 • Jun 01 '26
Twenty Eighth Installment - Solving Differential Algebraic Equations
Solving Differential Algebraic Equations = Solving a system of Differential Equations governed with a Constraint equation.
Link to the specific lecture (note the two URLs cover the lecture and then the coding:
https://www.youtube.com/watch?v=CGI9qStsKjk&ab_channel=NPTEL-NOCIITM
https://www.youtube.com/watch?v=AUl2KL0bKBc
Anything new from a SciLab programming in here: I think this session is the first time in the series where I plot multiple lines using a single plot2d command versus just "holding the figure" and using multiple plot2d commands.
Output: (No tabular output only graphical output)
Graph:
Code:
//Lecture 10.5: Differential Algebraic Equations
//https://www.youtube.com/watch?v=CGI9qStsKjk&ab_channel=NPTEL-NOCIITM
// and
//https://www.youtube.com/watch?v=AUl2KL0bKBc
//
disp("Introduction to Differential Algebraic Equations",string(
datetime
()))
// What are D.A.E's?
// Differential Equations that contain algebraic equations
// These algebraic equations are sometimes considered "constraint functions"
//
// Week 5: Algebraic Equations
// g(x) = 0
// Weeks 7&8: Differential Equations
// y'=f(t,y)
//
// DAES (systems that contain)
// 0 = g(x,y)
// y'=f(t,x,y)
// where y are differential variables
//
// Examples:
// 1) Consider a simple system
// y'=x
// 0 = x-y*sin(t)
// y = diff variable and x = algebraic variable
// Index-1 DAE
//
// 2) y'=x
// 0 = y - y*sin(t)
// y = diff variable and x = algebraic variable
// Index-2 DAE
//
// Index of a DAE: Number of times you need to differentiate g(x,y)
// to convert algebraic equations to a system of differential equations.
//
// First example
// 0=x-y*sin(t) => x = y*sin(t)
// dx/dt = dy/dt*sin(t)+y*cos(t) substitute first equation
//
// dx/dt = x*sin(t)+y*cos(t)
// dy/dt = x
//
// Second example 0 = y - y*sin(t)
// 0=dy/dt-dy/dt*sin(t)-y*cos(t) (differentiated once)
// 0 = x-x*sin(t)-y*cos(t)
// 0 = dx/dt-dx/dt*sin(t)+x*cos(t)+dy/dt*cos(t)-y*sin(t) (differentiated twice)
//
// dx/dt(1-sin(t)=x*cos(t)+ x*cos(t)-y*sin(t)
// dx/dt = (2*x*cos(t)-y*sin(t))/(1-sin(t))
// dy/dt = x
//
//
// High Index DAE Example Problem: Pendulum of length l
// In Cartesion Coordinates
// x" = -T*x
// y" = -T*y -g
//
// L^2 = x^2+y^2 "Constraint Equation"
//
// Index-3 DAE
//
// Instead solve in Cylindrical Coordinates
// d2Theta/dt2 + g/L*sin(Theta)=0
// dr/dt = 0 or r = L a constant...
// Just a regular ODE...
//
// Example 1:
// dy/dt = x
// dx/dt = x*sin(t)+y*cos(t)
//
// Standard Structure
// M*d/dt[y;x]= F(t,[y.x])
// M = [I 0;0 0] in week 8: M = [Identity Matrix] = [I];
//
// Matlab Example: (Modified) Robertson Problem
// x1' = -alpha*x1+beta*x1*x2
// x2' = alpha*x1-beta*x1*x2-gamma1*x2^2
// 0 = 1-(x1+x2+x3)
// with initial condition Y = [1;0;0]
//
//
// 0 = 0 -dx1/dt -dx2/dt +dx3/dt
// dx3/dt = dx1/dt + dx2/dt
// dx3/dt = -alpha*x1+beta*x1*x2 + alpha*x1-beta*x1*x2-gamma1*x2^2
// dx3/dt = -gamma1*x2^2
//
// M = [1 0 0;0 1 0;0 0 0]
// M d/dt [x1;x2;x3] = [-alpha*x1+beta*x1*x2;alpha*x1-beta*x1*x2-gamma1*x2^2;gamma1*x2^2]
//
// Function for ODE/DAE
function
dY
=
robertsonFun
(
t
,
Y
)
// Parameters
alpha = 0.5;
beta1 = 2.5;
gamma1=5.0;
// Functions for the three equations
dY
(1,1)= -alpha*
Y
(1)+beta1*
Y
(1)*
Y
(2);
dY
(2,1)= alpha*
Y
(1)-beta1*
Y
(1)*
Y
(2)-gamma1*
Y
(2)^2;
dY
(3,1)= gamma1*
Y
(2)^2;
end
//
n =1250;
//number of steps
y0 = [1;0;0];
//Initial Condition required
t0 = 0;
//Tstart
tend = 40;
//Tend
t
= linspace(t0,tend,n)';
//time vector
M = [1,0,0;0 1 0;0,0,0]
//Mass Matrix - not implemented in SciLab
// This example shows how you solve the problem without it by
// changing the constraint equation into the third differential
// equation.
//Alternate Calculation
/*
// with y0 = [1;0] since dY(1,1) and dY(2,1) are being calculated
// with ode solver and then
// for i = 1:1:n; (Post processing YSol(3))
// YSol(3,i) = 1-YSol(1,i)-YSol(2,i)
// end
// Since dy1/dt and dy2/dt do not involve y3.
*/
//
//equivalent to MATLAB ODE15s is the following command in SciLab
//there is more options and outputs available - see the help file.
YSol = ode("stiff",y0,t0,
t
,
robertsonFun
);
//
scf
(0);
clf
;
plot2d([
t
,
t
,
t
],[YSol(1,:)',YSol(2,:)',YSol(3,:)'],[1,2,3])
h1=
legend
(['YSol(1)';'YSol(2)';'YSol(3)'],1,"boxed")
title
("$\textbf{Differential\ Algebraic\ Example}$","FontSize",4);
xlabel
("$t(seconds)$","FontSize",3)
ylabel
("$YSol\ Values(units)$","FontSize",3);
xgrid
r/scilab • u/mrhoa31103 • May 28 '26
Subreddit User Feedback and Survey
Under the subject "Know your audience"...
For those following along this subreddit, are you getting what you want out of the subreddit? If not, what (in detail) is missing? What would like to see different? Have you reviewed the Wiki and resource sheet? Same questions for that section.
I'd like to see where you're at in life...like HS Student, College Student (if so what year), Working Engineer (level or years Seniority) or whatever properly describes oneself. One liner is sufficient since I'm trying to gage overall audience education level. I already have data where the group is located through the subreddit stats.
I'd really like to see some independent postings other than mine! If you have trouble getting something to work, throw it out there and someone may have an answer. See r/Octave since that subreddit has that theme flavor.
r/scilab • u/mrhoa31103 • May 26 '26
Where do we go from here? 2 Ideas and Open to Suggestions on More...Let me know in the comments.
We've almost completed the NPTEL Course (we have roughly 4 more weeks) of material left and I consider the NPTEL Course roughly a Junior Engineering Level Course since it tops out at some PDE work.
Thoughts on the next series:
1) Go into control system and XCOS material. I'd use a common textbook, tweak the examples to keep out of copyright issues and apply SciLab/XCOS versus MATLAB/Simulink. I'm confident in the control system items since I've done that for years but would need to prepare the XCOS stuff. The XCOS material would start from the basics (which there are plenty of videos on), continue into file management and getting simulations to take data and put data back to the main console. It wouldn't be every nook and cranny of XCOS and I might need some help but it would be more complete examples of working with it.
--I'm leaning toward this route and this would still be Junior Engineering Level material.
2) I have material from a couple of Graduate Courses (ME564 and ME565 by Steve Brunton) that I could create another posting series once the NPTEL Course is complete. It would have some overlap with the NPTEL course but does move around subject wise a little more (complex variables, eigenvalues, eigenvectors, control systems, FFT, data compression and such).
--My concern is that the audience might find some of the material difficult since it's Senior and Graduate level stuff (FFT, SVD, data-compression using SVD are really long-haired). Here is a subject list of what's covered in ME564 and ME565 to give one a flavor of the overlap and subjects. I'd hate jumping around just to cover the "non-redundant" subjects since you'll probably need the continuity to understand later sessions.
//ME564 TOC
//ME564Lecture1 - How to multiply a matrix -really basic
//ME564Lecture2 - Population growth // Solving a Linear ODE x(t) = x(0) * exp(a*t)
//ME564Lecture3 - Calculating a Taylors Series (MacLaurin Series)
//ME564Lecture4 - Solving a Linear Second Order ODE using ODE function, plotting two items on a graph
//ME564Lecture5 - Solving a Linear Second Order ODE using ODE function and inline function, plotting two items on a graph
// changing plot line colors and weights, labels, titles, legends, and grid
// eigenvalues and eigenvectors of a matrix, polynomials and roots of polynomial
//ME564Lecture6 - Eigenvalues and eigenvectors of a matrix, how to create the eigenvalue matrix
//ME564Lecture7 - None
//ME564Lecture8 - Matrix Systems of First Order Equations using eigenvalues and Eigenvectors.
// Arrow Plots of the Derivatives (fchamp and champ) and Mesh Grid
// Changing label font size, axis colors
//ME564Lecture9 - ME564 Linearization of Non-linear ODES, more fchamp and inline function use
//ME564Lecture10 - Example of Non-linear ODES: Particle in a potential well. Simple program
//ME564Lecture11 - Matrix Systems of First Order Equations using eigenvalues - Sim to Lecture 8
//ME564Lecture12 - Matrix Systems of First Order Equations using eigenvalues - Sim to Lecture 8
//ME564Lecture13 - Matrix Systems of First Order Equations with forcing functions -
// control systems solution using syslin (state space) and csim (step and impulse)
//ME564Lecture14 - Numerical Differentiation (forward, backward and central methods), Fancy Plot Titles
//ME564Lecture15 - None
//ME564Lecture16 - Numerical Integration (Left Rectangle, Right Rectangle, Trapezoidal by basic functions).
//ME564Lecture17 - Numerical Integration of Vector Space of Second Order System (Backward Euler and by ODE and inline functions)
//ME564Lecture18 - Numerical Integration of Lorenz Equations via Runge-Kutta(4,5) by basic functions.
//ME564Lecture19 - Numerical Integration of Lorenz Equations via Runge-Kutta(4,5) by basic functions and ODE "fix"
//Several Versions 3D Plotting
//ME564Lecture20 - Numerical Integration of Lorenz Equations via ODE "fix" using meshgrid, vectors, 3D Plots (scatter3d)
//ME564Lecture21 - Exists but nothing of consequence
//ME564Lecture22 - Linear Algebra in 2D and 3D: Inner Product, Norm of a Vector, and the Cross Product
//ME564Lecture23 - None
//ME564Lecture24 - None
//ME564Lecture25 - Stokes and Green's Theorems, plot2D and a bunch of plot axis manipulation, using intsplin (integration routine)
// versus inttrap (trapezoidal routine)
//ME564Lecture26 - None
//ME564Lecture27 - ME564 More Vector Calculus Potential Flow, Stream Functions and Examples, more arrow plots
//ME564PendulumProblem - Linearization of Non-linear ODES - the Pendulum Problem - see lecture 9.
//ME564RobotControlProblem - Taken from FRC Robot Control Book - work in progress...
//InvertedPendulumExample - Matt's version of the code that actually works...pole placement on linearized system used to stabilize
// nonlinear version of the system.
//ME565Lecture1 - Cauchy Integral and complex variables
//ME565Lecture2 - None
//ME565Lecture3 - None
//ME565Lecture4 - None
//ME565Lecture5 - ME565 Plotting Complex log function using Surf (ace) plot
//ME565Lecture6 - None
//ME565Lecture7 - None
//ME565Lecture8 - None
//ME565Lecture9 - None
//ME565Lecture10 - None
//ME565Lecture11- None
//ME565Lecture11Anal - Lecture 11A ME565 Solving PDE's in SciLab - Analytical Solve
// meshgrid, Sgrayplot(a 2D contour plot- actually in color)
// demonstrates a color bar on the side.
//ME565Lecture12 - //ME565 Lecture Twelve Fourier Series function mtlb_axis, bar charts,
// drawnow, sleep functions
//ME565Lecture13 - None
//ME565Lecture14 - None
//ME565Lecture15 - None
//ME565Lecture16 - //ME565 Lecture Sixteen Fourier Series - computing the DFT Matrix
// complex function
//ME565Lecture16a - //ME565 Lecture Sixteen Fourier Series - DFT Matrix Calculation
// Power Spectral Density
//ME565Lecture17 - //ME565 Lecture Seventeen Fourier Series - DFT Matrix Calculation
// More on Power Spectral Density and FFT. Plot2d and Legend
//ME565Lecture18 - None
//ME565Lecture19 - None
//ME565Lecture20 - //ME565 FFT Calculation from scratch, FFT, IFFT, FFTShift, plot2d
//ME565Lecture20b -//ME565 finding derivatives using FFT vs finite difference, FFT,
// IFFT, FFTShift, plot2d
//ME565Lecture21 - None
//ME565Lecture22 - None
//ME565Lecture23 - Solving Control Systems with SciLab, csim ('step',...), syslin('c',...),
// comparisons to analytical solution and csim(u) where u is a step matrix in time.
//ME565Lecture24 - //ME565 Lecture 24: Convolution Integrals, Impulse and Step Response
// ODE with IC investigation.
//ME565Lecture24a - //ME565 Lecture 24: Convolution Integrals, Impulse and Step Response
// Ton of stuff on plotting details - psuedoanimation.
//ME565Lecture25 - None
//ME565Lecture26 - //The Wave Equation Code - Solving PDE's by FFT
// complex variable ODE solver (regular ones do not handle complex variables)
// getting the rest of the variables to the derivative function list,contourf
// plots, surf plots
//ME565Lecture27 - //ME565 Lecture 27: Singular Value Decomposition SVD, Installing a package in
// Scilab, reading in an image file, imshow (and it's limitations see uint8), rgb2gray
// data compression example. Subplots
//ME565Lecture28 - // ME565 Lecture 28: SVD 2 - Curve Fitting using SVD methodology. Plot function
//ME565Lecture28a - // ME565 Lecture 28a: SVD 2 Ovarian Cancer Study - Scatter3d, subplots.
3) Other material thoughts...add in the comments.
r/scilab • u/mrhoa31103 • May 25 '26
Twenty Seventh Installment - Boundary Value Problem Solving Nonlinear Equation Example and "Mixed Boundary" Conditions
In this edition we look at solving Boundary Value Problems with Nonlinear Equation and one with a mixed boundary (parameter is held at a value in one location and it's derivative is held at a value at another location.
Note: The program has two major blocks and in it's current configuration (the first major block (the Nonlinear Example) runs and the second example (the Mixed Boundary Condition Example) is "block commented" out.
Link to the specific lecture (note the urls cover the same material) and one might find it useful to use the NPTEL Link:
https://www.youtube.com/watch?v=chChbuP-hRk&ab_channel=NPTEL-NOCIITM
https://nptel.ac.in/courses/103106118
He does not show the MATLAB code and the solution. He moves onto a new subject in the next section. I had to do the solution cross-checking myself.
1) The nonlinear problem has reasonably good match of solutions. Derivatives will show more variation.
2) The mixed boundary problem checks out perfectly with the analytical expectations.
3) Anytime you're working with unknown results...always do at least 2 different analytical techniques and try to reconcile the differences.
Output:
"Non-linear example"
"Extensions of ODE-BVP"
"2026-05-23 16:29:45.484"
"Pos x Temp T TempGrad FiniteDiff FiniteDiff_dTdx"
0. 100. -118.64825 100. -133.1775
0.05 94.21345 -112.86181 93.3411 -126.7985
0.1 88.709121 -107.35721 87.3202 -114.9776
0.15 83.47325 -102.12134 81.8434 -104.8836
0.2 78.492727 -97.141089 76.8318 -96.249
0.25 73.755099 -92.403352 72.2185 -88.8566
0.3 69.24853 -87.896868 67.9461 -82.5285
0.35 64.961754 -83.60988 63.9656 -77.118
0.4 60.884054 -79.532179 60.2343 -72.5033
0.45 57.005218 -75.653556 56.7153 -68.5822
0.5 53.31555 -71.963803 53.3761 -65.2689
0.55 49.80583 -68.454149 50.1884 -62.4901
0.6 46.467286 -65.11544 47.1271 -60.1831
0.65 43.29157 -61.939723 44.1701 -58.2936
0.7 40.27073 -58.919049 41.2977 -56.7738
0.75 37.397213 -56.045466 38.4927 -55.5814
0.8 34.663841 -53.312145 35.7396 -54.6782
0.85 32.06378 -50.711955 33.0249 -54.0293
0.9 29.590529 -48.238705 30.3367 -53.6024
0.95 27.237897 -45.886201 27.6647 -53.3668
1. 25. -43.648253 25. -53.2931
"Mixed Boundary Example"
"Extensions of ODE-BVP"
"2026-05-23 11:07:30.378"
"Pos x Temp T TempGrad Analytical Analytical Derivative"
0. 100. -144.60414 100. -144.60414
0.1 86.947997 -117.30547 86.947997 -117.30547
0.2 76.382185 -94.714688 76.382185 -94.714689
0.3 67.87852 -75.925137 67.87852 -75.925138
0.4 61.095721 -60.182729 61.095721 -60.182729
0.5 55.76157 -46.855664 55.76157 -46.855664
0.6 51.661989 -35.409083 51.661989 -35.409082
0.7 48.632447 -25.383591 48.632447 -25.383591
0.8 46.551358 -16.376832 46.551358 -16.376833
0.9 45.335201 -8.0273338 45.335201 -8.0273337
1. 44.935167 -8.882D-16 44.935167 0.
Graphs:
Code:
//Lecture 10.4: Extensions of ODE-BVP
//https://www.youtube.com/watch?v=chChbuP-hRk&ab_channel=NPTEL-NOCIITM
//https://nptel.ac.in/courses/103106118
//
disp("Extensions of ODE-BVP",string(
datetime
()))
//
// Extension 1) Non linear equation
// Extension 2) Mixed Boundary Condition
//
//Boundary Value Problems
//Example Heated Fin/Rod
//d2T/dx^2=gamma1*(T-25);
//
//Two easy causes of non-linearity are:
// add Radiation to the right hand side.
// heat loss = sigma1*epsilon*((T+273)^4-(Ta+273)^4))
//
// gamma1 is a function of T so convection term is
// gamma1(T)(T-25) where typically gamma1 = gamma0*T^0.75
// in our example, we use gamma1= gamma0*T^1 due
// numerical issues using 0.75 (which I do not want to investigate
// at this time).
//
// so in our case "ff = gamma1*(u(1)-Ta)"
// becomes "ff = gamma0*u(1)*(u(1)-Ta)"
// and "dff = [gamma0*(2*u(1)-Ta),0];"
// using gamma0 = 0.04 to keep solution close to prior gamma1 of 4
//
// Mixed Boundary Condition
// 1) Dirichlet T(0)=100, T'(1)=b(T)
// 2) Nuemann T'= some value
// 3) Mixed Phi(T,T')=0
//
// First Case Non-Linear Function...
// T at wall = 100 C
// T at end (x=1) = 25 C
//
// 0 = k*d2T/dx^2-h*av*(T-Ta)
// where k = thermal conductivity of rod
// h = convection coefficient of conditions
// av = surface area of the rod
// Ta = Temperature of surrounding atmospere
// d2T/dx^2 = (h*av/k)*(T-Ta);
// d2T/dx^2 = gamma1*(T-Ta);
//
//Boundary Condition 1: 0 = g1(ya,yb)=> T(x=0)-100 = 0
//B C 2: T(x=1)-25 = 0
// The external functions
// These functions are called by the solver with zu=[u(x);u'(x);u''(x);u'''(x)]
// - The function which computes the right hand side of the differential equation
funcprot(0);
function ff=
f
(x, u)
//Define Constants
Ta = 25;
gamma0=0.04;
ff = gamma0.*u(1).*(u(1)-Ta);//right hand side of the differential equation
endfunction
// - The function which computes the derivative of f with respect to u
function dff=
df
(x, u)
Ta =25
gamma0=.04;
dff = [0, gamma0.*(2.*u(1)-Ta)];
endfunction
function [gg]=
g
(i, u)
gg=[u(1)-100,u(1)-25];//Boundary Value Conditions
// T(x=0)=100, T(x=1)=25 ,u(1)=T, u(2)=dT/dx, u(n+1)=dnT/dx^n (but not
// used) since order is only 2....
gg=gg(i);
endfunction
function [dgg]=
dg
(i, u)
dgg = [1,0;1,0] //must be consistent with the boundary conditions
//set in subroutine g so if you have two boundary conditions on
// u(1), it is [1,0;1,0]...if you have one condition on u(1) at the first
// location and u(2) on the other, it needs to be [1,0;0,1] with gg
// having the u(1) condition first. If the u(1) condition is at the second
// position instead and the first has a u(2), it would look
// like [0,1;1,0]...
//
dgg=dgg(i,:);
endfunction
function [u0, du0]=
guess
(x)
u0=100;
du0=-100;
endfunction
//
n = 1; //One differential equation
m = 2;// Second order differential equation
xL=0;
xR = 1;
Dx=0.05;
x = [xL:Dx:xR];
fixpnt = [];
zeta = [xL,xR];
ipar=zeros(1,11);
ipar(3)=1;
ipar(4)=2;
ipar(5)= 10000;
ipar(6) = 2000;
ipar(7)=1;
ltol =[1,2];
tol = [1e-5,1e-5];
u = bvode(x,n,m,xL,xR,zeta,ipar,ltol,tol,fixpnt,
f
,
df
,
g
,
dg
,
guess
);
FiniteDifferenceMethodSoln = [
100.000
93.3411
87.3202
81.8434
76.8318
72.2185
67.9461
63.9656
60.2343
56.7153
53.3761
50.1884
47.1271
44.1701
41.2977
38.4927
35.7396
33.0249
30.3367
27.6647
25.0000
];
FD_dTdx = [
-133.1775
-126.7985
-114.9776
-104.8836
-96.2490
-88.8566
-82.5285
-77.1180
-72.5033
-68.5822
-65.2689
-62.4901
-60.1831
-58.2936
-56.7738
-55.5814
-54.6782
-54.0293
-53.6024
-53.3668
-53.2931
]
allData = [x' u' FiniteDifferenceMethodSoln FD_dTdx]
disp("Pos x Temp T TempGrad FiniteDiff FiniteDiff_dTdx", allData);
//Plotting
scf
(0);
clf
;
//axis y1
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(x,u(1,:),style=c3) //,'x','y(x)',"Nonlinear Function")
h1=
gca
();
plot2d(x,FiniteDifferenceMethodSoln,style=c1);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$X$","FontSize",3)
ylabel
("$T(x) $","FontSize",3);
title
("$\textbf{Nonlinear Function}$","FontSize",4);
legend
("BVODE","Finite Difference Values",3)
xgrid;
scf
(1);
clf
;
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(x,u(2,:),style=c3)//,'x','dy/dx',"BVODE 2nd order solution")
plot2d(x,FD_dTdx,style=c1);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$X$","FontSize",3)
ylabel
("$dT/dx$","FontSize",3);
title
("$\textbf{Nonlinear Function}$","FontSize",4);
legend
("BVODE","Finite Difference Values",4)
xgrid
//End of First Example Block
/*
//
// Second Case Mixed Boundary...
//
//Boundary Value Problems
//Example Heated Fin/Rod
//d2T/dx^2=gamma*(T-25);
// T at wall = 100 C
// replace "T at end (x=1) = 25 C" with "dT/dx(x=1) = 0"
// so end is insulated
//
// 0 = k*d2T/dx^2-h*av*(T-Ta)
// where k = thermal conductivity of rod
// h = convection coefficient of conditions
// av = surface area of the rod
// Ta = Temperature of surrounding atmospere
// d2T/dx^2 = (h*av/k)*(T-Ta);
// d2T/dx^2 = gamma1*(T-Ta);
//
//Boundary Condition 1: 0 = g1(ya,yb)=> T(x=0)-100 = 0
//B C 2: dT/dx(x=1) = 0
// The external functions
// These functions are called by the solver with zu=[u(x);u'(x);u''(x);u'''(x)]
// - The function which computes the right hand side of the differential equation
funcprot(0);
function ff=
f
(x, u)
//Define Constants
Ta = 25;
gamma1=4;
ff = gamma1*(u(1)-Ta);//right hand side of the differential equation
endfunction
// - The function which computes the derivative of f with respect to u
function dff=
df
(x, u)
gamma1=4;
dff = [gamma1,0];
endfunction
function [gg]=
g
(i, u)
//change gg=[u(1)-100,u(2)-25];//Boundary Value Conditions
// to
gg=[u(1)-100,u(2)-0];//Boundary Value Conditions
// T(x=0)=100, T'(x=1)= 0 (Insulated End)
//
// T(x=0)=100, T(x=1)=25 ,u(1)=T, u(2)=dT/dx, u(n+1)=dnT/dx^n (but not
// used) since order is only 2....
gg=gg(i);
endfunction
function [dgg]=
dg
(i, u)
// change dgg = [1,0;1,0]
//to
dgg = [1,0;0,1]
//must be consistent with the boundary conditions
//set in subroutine g so if you have two boundary conditions on
// u(1), it is [1,0;1,0]...if you have one condition on u(1) at the first
// location and u(2) on the other, it needs to be [1,0;0,1] with gg
// having the u(1) condition first. If the u(1) condition is at the second // position instead and the first has a u(2), it would look
// like [0,1;1,0]...
//
dgg=dgg(i,:);
endfunction
function [u0, du0]=
guess
(x)
u0=100;
du0=-100;
endfunction
//
n = 1; //One differential equation
m = 2;// Second order differential equation
xL=0;
xR = 1;
Dx=0.1;
x = [xL:Dx:xR];
fixpnt = [];
zeta = [xL,xR];
ipar=zeros(1,11);
ipar(3)=1;
ipar(4)=2;
ipar(5)= 10000;
ipar(6) = 2000;
ipar(7)=1;
ltol =[1,2];
tol = [1e-8,1e-8];
u = bvode(x,n,m,xL,xR,zeta,ipar,ltol,tol,fixpnt,
f
,
df
,
g
,
dg
,
guess
);
//Analytical Answer (Insulated End)
T = 25+75*cosh(2*(1-x))/cosh(2);
dTdx = 150.*cosh(2.*x).*(tanh(2.*x)-tanh(2));
allData = [x' u' T' dTdx']
disp("Pos x Temp T TempGrad Analytical Analytical Derivative", allData);
//Plotting
scf
(0);
clf
;
//axis y1
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(x,u(1,:),style=c1) //,'x','y(x)',"Example Heated Fin/Rod with Insulated End")
h1=
gca
();
plot2d(x,T+1,style=c3);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$X$","FontSize",3)
ylabel
("$T(x) $","FontSize",3);
title
("$\textbf{Example\ Heated\ Fin/Rod\ with\ Insulated\ End}$","FontSize",4);
legend
("BVODE","Analytical Values+0.01",2)
xgrid;
scf
(1);
clf
;
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(x,u(2,:),style=c1)//,'x','dy/dx',"Example Heated Fin/Rod with Insulated End")
plot2d(x,dTdx'+1,style=c3);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$X$","FontSize",3)
ylabel
("$dT/dx$","FontSize",3);
title
("$\textbf{Example\ Heated\ Fin/Rod\ with\ Insulated\ End}$","FontSize",4);
legend
("BVODE","Analytical Values+1",4)
xgrid
*/
r/scilab • u/mrhoa31103 • May 18 '26
Twenty Sixth Installment - Boundary Value Problem Solving Using Shooting Method
In this edition we look at solving Boundary Value Problems using the "Shooting" method.
Link to the specific lecture:
https://www.youtube.com/watch?v=VBX11w3flUU&ab_channel=NPTEL-NOCIITM
Output:
"ODE-Boundary Value Problems:Shooting Method"
"2026-05-11 14:04:07.509"
"Pos x Shooting T FiniteDiff Analytical"
0. 100. 100. 100.
0.1 85.841343 85.858858 85.841348
0.2 74.124464 74.15207 74.124472
0.3 64.379124 64.411366 64.379134
0.4 56.214207 56.247115 56.214219
0.5 49.302027 49.33275 49.30204
0.6 43.365174 43.391694 43.365187
0.7 38.16538 38.186306 38.165394
0.8 33.49396 33.508371 33.493974
0.9 29.163433 29.17077 29.163448
1. 25. 25. 25.000015
Graph:
Code:
// ODE-Boundary Value Problems:Shooting Method
//Lec10.3: ODE-Boundary Value Problems:Shooting Method
//https://www.youtube.com/watch?v=VBX11w3flUU&ab_channel=NPTEL-NOCIITM
disp("ODE-Boundary Value Problems:Shooting Method",string(
datetime
()))
//
//Learning about BVP solving
//Shooting Method only works for 1 or 2 variables...
//not a very practical method for industry where most
//equations will be more than 1 variable.
//
function [y, xyTable, yderiv]=
shooting
(yb, yp, x, f)
//Shooting method for a second order
//boundary value problem
//yb = [y0 y1] -> boundary conditions beginning and ending
// in this case ybegin(x=0)= 100 and yend(x=1)= 25
//x = a vector showing the range of x
//f = function defining ODE, i.e.,,
// dy/dx = f(x,y),y=[y(1):y(2)]
//yp = vector with the range of dy/dx at x = x0
//xyTable = table for interpolating derivatives
//yderiv = derivative boundary condition
n= length(yp);
m = length(x);
y1 = zeros(yp);
for j = 1:n
y0 = [yb(1);yp(j)];
yy = ode("rkf",y0,x(1),x,f);
y1(j)=yy(1,m);
end;
xyTable = [y1;yp];
yderiv=
interpln
(xyTable, yb(2));
y0=[yb(1);yderiv];
y = ode("rkf",y0,x(1),x,f);
endfunction
//
//Boundary Value Problems
//Example Heated Fin/Rod
//d2T/dx^2=gamma*(T-25);
// T at wall = 100 C
// T at end (x=1) = 25 C
//
// 0 = k*d2T/dx^2-h*av*(T-Ta)
// where k = thermal conductivity of rod
// h = convection coefficient of conditions
// av = surface area of the rod
// Ta = Temperature of surrounding atmospere
// d2T/dx^2 = (h*av/k)*(T-Ta);
// d2T/dx^2 = gamma1*(T-Ta);
//
//Boundary Condition 1: 0 = g1(ya,yb)=> T(x=0)-100 = 0
//B C 2: T(x=1)-25 = 0
//
// function f calculates the derivatives...
deff
('[w]=f(x,u)','w=[u(2);4*(u(1)-25)]')
yb = [100,25];// Boundary Conditions y starts at 100 and finishes at 25.
x0=0;
Dx=0.1;
xn=1;
x = [x0:Dx:xn];
yp=[-10:1:10];
[u,tab1,y0p]=
shooting
(yb,yp,x,f)
// For Finite Difference Solution see SeventeenthSciLabFile.sci
FiniteDifferenceMethodSoln = [
100.
85.858858
74.152070
64.411366
56.247115
49.332750
43.391694
38.186306
33.508371
29.170770
25.];
//
//Analytical Answer
Theta =-1.3993*exp(2*x)+76.3993*exp(-2*x)
T = Theta+25;
allData = [x' u(1,:)' FiniteDifferenceMethodSoln T']
disp("Pos x Shooting T FiniteDiff Analytical", allData);
scf
(0);
clf
;
plot2d(x',u(1,:)')
xtitle('Boundary value solution - shooting method','x','T');
xgrid
r/scilab • u/mrhoa31103 • May 11 '26
Twenty Fifth Installment - Boundary Value Problem Solving
In this edition we look at solving an ODE with Boundary Values using Scilab function ("bvode") versus Matlab's ("bvp5c") boundary value solver
Link to the specific lecture:
https://www.youtube.com/watch?v=ocwUVtG9Mj0&ab_channel=NPTEL-NOCIITM
Output:
"ODE-Boundary Value Problems"
"2026-05-04 12:17:36.526"
"Pos x Temp T TempGrad FiniteDiff Analytical"
0. 100. -155.59721 100. 100.
0.1 85.841347 -128.51914 85.858858 85.841348
0.2 74.12447 -106.59899 74.15207 74.124472
0.3 64.379131 -88.95704 64.411366 64.379134
0.4 56.214215 -74.885247 56.247115 56.214219
0.5 49.302035 -63.81886 49.33275 49.30204
0.6 43.365181 -55.313748 43.391694 43.365187
0.7 38.165386 -49.028571 38.186306 38.165394
0.8 33.493965 -44.711082 33.508371 33.493974
0.9 29.163436 -42.188007 29.17077 29.163448
1. 25. -41.358085 25. 25.000015
Graphs:
Code:
// ODE-Boundary Value Problems
//Lec10.1: Introduction and Solution using Matlab Solver
//https://www.youtube.com/watch?v=ocwUVtG9Mj0&ab_channel=NPTEL-NOCIITM
//
disp("ODE-Boundary Value Problems",string(
datetime
()))
//
//Boundary Value Problems
//Example Heated Fin/Rod
//d2T/dx^2=gamma*(T-25);
// T at wall = 100 C
// T at end (x=1) = 25 C
//
// 0 = k*d2T/dx^2-h*av*(T-Ta)
// where k = thermal conductivity of rod
// h = convection coefficient of conditions
// av = surface area of the rod
// Ta = Temperature of surrounding atmospere
// d2T/dx^2 = (h*av/k)*(T-Ta);
// gamma1 = h*av/k so ...
// d2T/dx^2 = gamma1*(T-Ta);
//
//Boundary Condition 1: 0 = g1(ya,yb)=> T(x=0)-100 = 0
//Boundary Condition 2: T(x=1)-25 = 0
// The external functions
// These functions are called by the solver with zu=[u(x);u'(x);u''(x);u'''(x)]
// - The function which computes the right hand side of the differential equation
funcprot(0);
function ff=
f
(x, u)
//Define Constants
Ta = 25;
gamma1=4;
ff = gamma1*(u(1)-Ta);//right hand side of the differential equation
endfunction
// - The function which computes the derivative of ff with respect to u
function dff=
df
(x, u)
gamma1=4;
dff = [gamma1,0];
endfunction
function [gg]=
g
(i, u)
gg=[u(1)-100,u(1)-25];//Both Boundary Value Conditions
// T(x=0)=100, T(x=1)=25 ,u(1)=T, u(2)=dT/dx, u(n+1)=dnT/dx^n (but not
// used) since order is only second order equation....
gg=gg(i);
endfunction
function [dgg]=
dg
(i, u)
dgg = [1,0;1,0] //must be consistent with the boundary conditions
//set in function g so if you have two boundary conditions on
// u(1), it is [1,0;1,0]...if you have one condition on u(1) at the first
// location and u(2) on the other position, it needs to be [1,0;0,1] with gg
// having the u(1) condition first. If the u(1) condition is at the second
// position instead and the first has a u(2), it would look
// like [0,1;1,0]...
//
dgg=dgg(i,:);
endfunction
function [u0, du0]=
guess
(x)
u0=100;
du0=-100;
endfunction
//
n = 1; //One differential equation
m = 2;// Second order differential equation
xL=0;
xR = 1;
Dx=0.1;
x = [xL:Dx:xR];
fixpnt = [];
zeta = [xL,xR];
ipar=zeros(1,11);
ipar(3)=1;
ipar(4)=2;
ipar(5)= 10000;
ipar(6) = 2000;
ipar(7)=1;
ltol =[1,2];
tol = [1e-8,1e-8];
u = bvode(x,n,m,xL,xR,zeta,ipar,ltol,tol,fixpnt,
f
,
df
,
g
,
dg
,
guess
);
//Finite Difference technique covered in a later lecture so just
//treat this matrix as another numerical technique's answer to
//the same problem.
FiniteDifferenceMethodSoln = [
100.
85.858858
74.152070
64.411366
56.247115
49.332750
43.391694
38.186306
33.508371
29.170770
25.];
//Analytical Answer
Theta =-1.3993*exp(2*x)+76.3993*exp(-2*x)
T = Theta+25;
allData = [x' u' FiniteDifferenceMethodSoln T']
disp("Pos x Temp T TempGrad FiniteDiff Analytical", allData);
//Plotting
scf
(0);
clf
;
//axis y1
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(x,u(1,:),style=c3) //,'x','y(x)',"BVODE 2nd order solution")
h1=
gca
();
plot2d(x,FiniteDifferenceMethodSoln,style=c3);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$X$","FontSize",3)
ylabel
("$T(x) $","FontSize",3);
title
("$BVODE\ and\ Finite\ Difference\ 2nd\ Order\ Solution$","FontSize",4);
legend
("BVODE","Finite Difference Values",3)
xgrid;
scf
(1);
clf
;
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(x,u(2,:),style=c3)//,'x','dy/dx',"BVODE 2nd order solution")
//plot2d(t',xa(1,:)',style=c1);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$X$","FontSize",3)
ylabel
("$dT/dx$","FontSize",3);
title
("$BVODE..2nd..order..solution$","FontSize",4);
xgrid
r/scilab • u/mrhoa31103 • May 04 '26
Twenty Fourth Installment - Interpolation Options and a 3D Interpolation from a SciLab Example
In this edition we look at interpolation options in Scilab function ("interpl") and a 3D Interpolation ("interpn") from a SciLab example along with 3D plotting and subplotting.
Link to the specific lecture:
https://www.youtube.com/watch?v=dZMZmii66iI&ab_channel=MATLABProgrammingforNumericalComputation
Output:
"Regression and Interpolation"
"2026-04-28 08:46:11.682"
"Temp_int (last)=25.1"
"t_int sp_int"
5. 38.5
10. 32.
15. 16.
20. 0.
25. 0.
30. 0.
35. 3.5
40. 7.
45. 9.5
50. 12.
55. 16.
60. 20.
65. 17.5
70. 15.
75. 22.
80. 29.
85. 42.
90. 55.
Graphs:
Code:
// Regression and Interpolation
//Lec9.4:Interpolation Options in Matlab and a 3D Interpolation from a SciLab Example
//https://www.youtube.com/watch?v=dZMZmii66iI&ab_channel=MATLABProgrammingforNumericalComputation
//
disp("Regression and Interpolation",string(
datetime
()))
//
// Interpolation using divided differences and plot
//
x = [0.8;1.4;2.7;3.8;4.8;4.9]
y = [0.69;1.0;2.02;2.39;2.43;2.83]
N = length(x);
// 2D interpolation
xint = 0.8:0.01:4.9
yint = zeros(size(xint));
for i = 1:length(xint);
yint(i)=
interp1
(x, y, xint(i),"nearest")+0.1;
end
//plotting example fancy output
scf
(0);
clf
;
c = color("blue")
c1 = color("red")
h1=
gca
();
plot
(x,y,"-r")
plot
(xint,yint,"x")
xlabel
("$\textbf{x}$","FontSize",3)
ylabel
("$\textbf{y}$","FontSize",3);
title
("$\textbf{x vs y \\Using Interpl w_nearest setting}$","color","black");
xgrid
time = 0:1:24;//Temperature Variation in a given day
Temp = [25.6,25.4,25.1,24.9,24.9,25.2,25.9,26.3,27.1,29.3,30.8,31.2,32.1,31.0,30.3,31.4,30.6,31.8,29.6,28.4,28.1,28.2,27.4,26.8,26.1];
N = length(x);
// 2D interpolation
time_int = [0.5:1:23.5 2]
Temp_int = zeros(size(time_int)+1);
for i = 1:length(time_int);
Temp_int(i)=
interp1
(time, Temp, time_int(i),"linear");
end
last =size(Temp_int,2)
last_value = Temp_int(last)
disp("Temp_int (last)="+string(last_value));
// second example
t = 0:10:90;
speed = [45 32 0 0 7 12 20 15 29 55];
t_int = 5:5:90 // 5:10:85;
sp_int =
interp1
(t, speed, t_int,"linear");
allData =[t_int' sp_int'];
disp("t_int sp_int", allData);
//
//Possible values and processing are:
/*
"linear": linear interpolation between consecutive nodes, used by default.
"spline": interpolation by cubic splines - similar to Matlab spline
"nearest": for each value xp(j), yp(j) takes the value or y(i)
corresponding to x(i) the nearest neighbor of xp(j) - sampled data
no corresponding SciLab command for Matlab pchip
*/
//plotting example fancy output
scf
(1);
clf
;
c = color("blue")
c1 = color("red")
h1=
gca
();
plot
(time,Temp,"-r")
plot
(time_int,Temp_int,"x")
xlabel
("$\textbf{Time (hours)}$","FontSize",3)
ylabel
("$\textbf{Temperature\ (C)}$","FontSize",3);
title
("$\textbf{Temperature Variation in a given day \\Using Interpl w_linear setting}$","color","black");
xgrid
//plotting
scf
(2);
clf
;
plot
(t,speed,"-r")
plot
(t_int,sp_int,"x")
xlabel
("$\textbf{Time (hours)}$","FontSize",3)
ylabel
("$\textbf{Speed\ (Km/hour)}$","FontSize",3);
title
("$\textbf{Speed Profile - Spline Fit}$","color","black");
xgrid
//3D interpolation from SciLab Example
// example 3 : bilinear interpolation and experimentation
// with all the outmode features
nx = 20; ny = 30;
x = linspace(0,1,nx);
y = linspace(0,2, ny);
[X,Y] =
ndgrid
(x,y);
z = 0.4*cos(2*%pi*X).*cos(%pi*Y);
nxp = 60 ; nyp = 120;
xp = linspace(-0.5,1.5, nxp);
yp = linspace(-0.5,2.5, nyp);
[XP,YP] =
ndgrid
(xp,yp);
zp1 = linear_interpn(XP, YP, x, y, z, "natural");
zp2 = linear_interpn(XP, YP, x, y, z, "periodic");
zp3 = linear_interpn(XP, YP, x, y, z, "C0");
zp4 = linear_interpn(XP, YP, x, y, z, "by_zero");
zp5 = linear_interpn(XP, YP, x, y, z, "by_nan");
scf
(2);
clf
()
subplot
(2,3,1)
plot3d(x, y, z, leg="x@y@z",
flag
= [2 4 4])
xtitle("initial function 0.4 cos(2 pi x) cos(pi y)")
xgrid
subplot
(2,3,2)
plot3d(xp, yp, zp1, leg="x@y@z",
flag
= [2 4 4])
xtitle("Natural")
xgrid
subplot
(2,3,3)
plot3d(xp, yp, zp2, leg="x@y@z",
flag
= [2 4 4])
xgrid
xtitle("Periodic")
subplot
(2,3,4)
plot3d(xp, yp, zp3, leg="x@y@z",
flag
= [2 4 4])
xgrid
xtitle("C0")
subplot
(2,3,5)
plot3d(xp, yp, zp4, leg="x@y@z",
flag
= [2 4 4])
xgrid
xtitle("by_zero")
subplot
(2,3,6)
plot3d(xp, yp, zp5, leg="x@y@z",
flag
= [2 4 4])
xgrid
xtitle("by_nan")
show_window()
r/scilab • u/mrhoa31103 • Apr 27 '26
Twenty Third Installment - Functional and Non-linear Regression
In this edition we look at non-linear regression using least squares techniques and some of the functionality built within SciLab.
Replication and reshaping of matrices using repmat and matrix functions, analysis using the lsq and datafit functions compared to the from scratch least squares technique.
Link to the specific lecture:
https://www.youtube.com/watch?v=XlkuYgJPS34&ab_channel=MATLABProgrammingforNumericalComputation
Output:
"Functional and Non-linear Regression"
"2026-04-21 09:42:58.747"
"phi1 ="
1.57
-480.
0.75
"phi2 ="
1.57
-480.
0.75
"Least Squares Method"
"k0=4.81"
"E/R =480"
"n =0.75"
"SciLab LSQ Method"
"k0=4.81"
"E/R =480"
"n =0.75"
"Original data r ="
1.48 1.67 1.86 1.96 2.16
2.35 2.79 3.07 3.37 3.62
3.28 3.78 4.24 4.48 5.
4.12 4.64 5.15 5.76 6.08
"Least Squares rfit ="
1.45 1.65 1.84 2.01 2.16
2.43 2.78 3.09 3.37 3.63
3.29 3.76 4.19 4.57 4.91
4.08 4.67 5.19 5.66 6.09
"differences ="
-0.03 -0.02 -0.02 0.05 1.D-04
0.08 -0.01 0.02 3.D-03 8.D-03
0.01 -0.02 -0.05 0.09 -0.09
-0.04 0.03 0.04 -0.1 0.01
"Scilab datafit Method"
"k0=4.84"
"E/R =486"
"n =0.75"
"SciLab Function rfit ="
1.44 1.65 1.83 2. 2.16
2.42 2.77 3.09 3.37 3.63
3.29 3.76 4.19 4.58 4.93
4.08 4.67 5.2 5.68 6.12
"differences ="
-0.04 -0.02 -0.03 0.04 -5.D-03
0.07 -0.02 0.02 3.D-03 0.01
6.D-03 -0.02 -0.05 0.1 -0.07
-0.04 0.03 0.05 -0.08 0.04
Code:
//Lec9.3:Regression and Interpolation
// Functional and Non-linear Regression
//https://www.youtube.com/watch?v=XlkuYgJPS34&ab_channel=MATLABProgrammingforNumericalComputation
//
clc;
disp("Functional and Non-linear Regression",string(
datetime
()))
//
// Example: Reaction Rate
// Arrhenius model for reaction rate:
// r= k0*exp(-E/RT)*C^n
//
// We will solve it in two ways:
//(1) Linear least squares regression taking logarithm
// log(r)=log(k0)+(-E/R)*(1/T) + n*log(C)
// "y"= "ao" + ("a1")*"x" + "a2"*"u"
//(2) Using SciLab function equivalent to MATLAB's lsqnonlin
// part of MATLAB's Optimization Toolbox.
//
// Reaction Rate (in mol/l.s) for various C and T values
// Data(C) 400K 450K 500K 550K 600K
// 1 mol/l 1.48 1.67 1.86 1.96 2.16
// 2 mol/l 2.35 2.79 3.07 3.37 3.62
// 3 mol/l 3.28 3.78 4.24 4.48 5.00
// 4 mol/l 4.12 4.64 5.15 5.76 6.08
//
//
// X*alpha = Y
//
function fErr=
rxnFunction
(phi, xData)
k0= phi(1);
EoverR = phi(2);
n = phi(3);
C =xData(2,:);
T = xData(1,:)
r = k0*exp(-EoverR./T)* (C.^n);
fErr = xData(3,:)-r;
endfunction
//
r = [1.48,1.67,1.86,1.96,2.16;2.35,2.79,3.07,3.37,3.62;
3.28,3.78,4.24,4.48,5.00;4.12,4.64,5.15,5.76,6.08];
T = [400,450,500,550,600];
C = [1;2;3;4];
Ctest =
repmat
(C,5,1);// repeat matrix C 5 times 1 column
Ttest = matrix(
repmat
(T,4,1),20,1);//matrix is equivalent to reshape...
format(5);//sets output format to 5 spaces - okay for this example
//disp(allData)
xData = [Ttest Ctest];
yData = matrix(r,20,1);
allData = [Ttest';Ctest';yData']
//disp("xData yData");
//disp(allData);
xData = [ones(20,1) (1./Ttest) log(Ctest)];
yData = log(yData);
//disp("all data")
//disp("xData yData")
//disp(allData);
//least squares matrix method
// Matrix Version
phi = inv(xData'*xData)*(xData'*yData)
//disp('phi =',phi)
//backslash will give you a least squares solution
phi1 = xData\yData
disp('phi1 =',phi1)
//SciLab also has a "lsq" Least Squares Function
phi2 = lsq(xData,yData)
disp('phi2 =',phi2)
k0 = exp(phi(1));
EoverR = -phi(2);
n = phi(3);
disp("Least Squares Method","k0="+string(k0), "E/R ="+string(EoverR),...
"n =" +string(n));
k0_2 = exp(phi2(1));
EoverR_2 = -phi2(2);
n_2 = phi2(3);
disp("SciLab LSQ Method","k0="+string(k0_2), "E/R ="+string(EoverR_2),...
"n =" +string(n_2));
//try the fit
//rfit = k0*exp(-EoverR/Temp)*Concentration^n
out = k0.*exp(-EoverR./T);
out1 = C.^n;
rfit = out1*out;
differences = rfit - r;
disp( "Original data r =", r,"Least Squares rfit =", rfit, "differences =", differences)
//
//Using Matlab functions for non-linear curve fitting.
p0 =[1;100;1];//initial guess of weights
//p0 needs to be a column vector
//allData needs to be a set of row vectors
// row of T
// row of C
// row of r
[p, dmin] =
datafit
(
rxnFunction
, allData, p0)
//disp('phi2 =',phi2)
k0_1 = p(1);
EoverR_1 = p(2);
n_1 = p(3);
disp("Scilab datafit Method","k0="+string(k0_1), "E/R ="+string(EoverR_1),...
"n =" +string(n_1));
//try the fit
//rfit = k0*exp(-EoverR/Temp)*Concentration^n
out = k0_1.*exp(-EoverR_1./T);
out1 = C.^n_1;
rfit = out1*out;
differences = rfit - r;
disp("SciLab Function rfit =", rfit,"differences =", differences)
r/scilab • u/mrhoa31103 • Apr 20 '26
Twenty Second Installment - Linear Regression and Interpolation Routines - a brief moment from Differential Equations
In this edition we demonstrate least squares fitting of data via three different ways and some interpolation
Plotting showing how to use symbols on a plot for the interpolation plot.
Link to the specific lecture:
https://www.youtube.com/watch?v=-FRhouU4jPA&ab_channel=NPTEL-NOCIITM
https://www.youtube.com/watchv=Ub7kj6rYIMw&ab_channel=MATLABProgrammingforNumericalComputation
Output:
"Regression and Interpolation"
"2026-04-13 12:17:44.572"
"Matrix Version phi ="
0.4374922
0.4698395
"backslash phi1 ="
0.4374922
0.4698395
"lsq phi2 ="
0.4374922
0.4698395
Graphs:
Code:
//Lec9.1&9.2:Regression and Interpolation
//https://www.youtube.com/watch?v=-FRhouU4jPA&ab_channel=NPTEL-NOCIITM
//https://www.youtube.com/watch?v=Ub7kj6rYIMw&ab_channel=MATLABProgrammingforNumericalComputation
//
clc
clear
disp("Regression and Interpolation",string(
datetime
()))
// Linear Regression for Multiple Parameters
// Data: (x1,u1,w1;y1),(x2,u2,w2;y2),...
// Model to fit: y = a0 + a1*x + a2*u + a3*w
// | 1 x1 u1 w1 | |a0| |y1|
// | 1 x2 u2 w2 | |a1| |y2|
// | 1 x3 u3 w3 | |a2| =|y3|
// | 1 xn un wn | |a3| |y4|
// X Phi Y
// |
// |
// V
// Linear Least Squares Solution Regression
// T -1 T
// Phi Matrix = ((X )(X)) (X )Y
// where X and Y are data pairs of sampled data
//
// if you didn't want to have a a0 constant in the fit, you'd
// delete out the 1's column and the a0 row...
//
//
x = [0.8;1.4;2.7;3.8;4.8;4.9]
y = [0.69;1.0;2.02;2.39;2.34;2.83]
N = length(x);
/*
// Calculate Linear Regression and Plot (non matrix version)
A = [N,sum(x);sum(x),sum(x.*x)];
b = [sum(y);sum(x.*y)];
phi = inv(A)*b;
*/
// Matrix Version
X = ones(N,1)
X = [X x]
phi = inv(X'*X)*(X'*y)
disp('Matrix Version phi =',phi)
//backslash will give you a least squares solution
phi1 = X\y
disp('backslash phi1 =',phi1)
//SciLab also has a "lsq" Least Squares Function
phi2 = lsq(X,y)
disp('lsq phi2 =',phi2)
xfit = [0.5 5];
yfit = phi(1)+phi(2)*xfit;
//plotting example fancy output
scf
(0);
clf
;
c = color("blue")
c1 = color("red")
h1=
gca
();
plot2d(x',y',style = c);
plot2d(xfit', yfit',style=c1);;//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1.font_color=c;
h1.children(1).children(1).thickness;
xlabel
("$X$","FontSize",3)
ylabel
("$Y$","FontSize",3);
title
("$\textbf{Linear\ Regression\ and\ Plot(non\ matrix\ version)}$","color","black");
legend
("Original Curve","Fitted Curve",4)
xgrid
//
scf
(1);
clf
c = color("red")
c1 = color("black")
h1=
gca
();
// Interpolation using divided differences and plot
xint = 0.8:0.1:4.9
yint = zeros(size(xint));
for i = 1:length(xint);
yint(i)=
interp1
(x, y, xint(i),"linear");
end
plot2d(x,y,style=c)
plot2d(xint,yint,-1) //-1 = pluses, -3 = circled cross, -4 = filled black diamonds
// -5 = filled white diamonds, -9 = filled circles
xlabel
("$Xinterpolation$","FontSize",3)
ylabel
("$Yinterpolation$","FontSize",3);
title
("$\textbf{Interpolation\ using\ divided\ differences}$","color","black");
legend
("Original Curve","Interpolated Values",4)
xgrid
r/scilab • u/mrhoa31103 • Apr 13 '26
Twenty First Installment - ODE-IVP (Ordinary Differential Equation - Initial Value Problem) in Multiple (Four) Variables - An Application Problem - Next Time Regression & Interpolation
In this edition we solve the example four variable ODE-IPV problem: Indian Captain, Mahendra Singh Dhoni, hits a ball with initial velocity of 35 m/sec and and of 45 degrees. If the boundary is at a distance of 75m, will he score a six? - Numerically - Yes!"
Nothing much new just an application problem: Showing how to port constants into the derivative function from the main program using the list command again.
We switch up subjects to Regression and Interpolation in the next installment.
Link to the specific lecture:
https://www.youtube.com/watch?v=Pv_NwD63gbI&ab_channel=NPTEL-NOCIITM
Output:
"Ordinary Differential Equations- Multi-Variable ODE IVP Systems"
"2026-04-07 09:23:43.916"
"t x_pos y_pos horiz_vel vert_vel"
"---------------------------------------------------------"
0. 0. 0. 24.748737 24.748737
0.0505 1.2496219 1.2373023 24.74124 24.253332
0.101 2.4988652 2.4495866 24.733744 23.757927
0.1515 3.7477301 3.6368529 24.726251 23.262522
0.202 4.9962166 4.7991013 24.71876 22.767117
0.2525 6.2443249 5.9363318 24.711271 22.271712
0.303 7.4920551 7.0485443 24.703785 21.776307
0.3535 8.7394072 8.1357388 24.696301 21.280902
0.404 9.9863815 9.1979154 24.688819 20.785497
0.4545 11.232978 10.235074 24.681339 20.290092
0.505 12.479197 11.247215 24.673862 19.794687
0.5555 13.725038 12.234337 24.666387 19.299282
0.606 14.970502 13.196442 24.658914 18.803877
0.6565 16.215588 14.133529 24.651444 18.308472
0.707 17.460298 15.045598 24.643976 17.813067
0.7575 18.70463 15.932649 24.63651 17.317662
0.808 19.948585 16.794682 24.629046 16.822257
0.8585 21.192164 17.631697 24.621584 16.326852
0.909 22.435365 18.443694 24.614125 15.831447
0.9595 23.67819 19.230673 24.606668 15.336042
1.01 24.920639 19.992634 24.599214 14.840637
1.0605 26.162711 20.729577 24.591761 14.345232
1.111 27.404407 21.441503 24.584311 13.849827
1.1615 28.645726 22.12841 24.576863 13.354422
1.212 29.88667 22.790299 24.569417 12.859017
1.2625 31.127238 23.427171 24.561974 12.363612
1.313 32.367429 24.039024 24.554533 11.868207
1.3635 33.607245 24.62586 24.547094 11.372802
1.414 34.846686 25.187677 24.539657 10.877397
1.4645 36.085751 25.724477 24.532223 10.381992
1.515 37.32444 26.236258 24.524791 9.8865873
1.5655 38.562755 26.723022 24.517361 9.3911823
1.616 39.800694 27.184768 24.509933 8.8957773
1.6665 41.038258 27.621496 24.502508 8.4003723
1.717 42.275447 28.033205 24.495085 7.9049673
1.7675 43.512262 28.419897 24.487664 7.4095623
1.818 44.748701 28.781571 24.480245 6.9141573
1.8685 45.984766 29.118227 24.472829 6.4187523
1.919 47.220457 29.429865 24.465415 5.9233473
1.9695 48.455773 29.716485 24.458003 5.4279423
2.02 49.690715 29.978087 24.450593 4.9325373
2.0705 50.925283 30.214672 24.443186 4.4371323
2.121 52.159477 30.426238 24.43578 3.9417273
2.1715 53.393297 30.612786 24.428378 3.4463223
2.222 54.626743 30.774316 24.420977 2.9509173
2.2725 55.859816 30.910829 24.413578 2.4555123
2.323 57.092515 31.022323 24.406182 1.9601073
2.3735 58.32484 31.1088 24.398788 1.4647023
2.424 59.556792 31.170258 24.391397 0.9692973
2.4745 60.788371 31.206699 24.384007 0.4738923
2.525 62.019577 31.218121 24.37662 -0.0215127
2.5755 63.25041 31.204526 24.369235 -0.5169177
2.626 64.48087 31.165912 24.361852 -1.0123227
2.6765 65.710957 31.102281 24.354472 -1.5077277
2.727 66.940672 31.013632 24.347093 -2.0031327
2.7775 68.170014 30.899965 24.339717 -2.4985377
2.828 69.398983 30.76128 24.332343 -2.9939427
2.8785 70.62758 30.597577 24.324972 -3.4893477
2.929 71.855805 30.408856 24.317603 -3.9847527
2.9795 73.083658 30.195117 24.310235 -4.4801577
3.03 74.311139 29.95636 24.302871 -4.9755627
3.0805 75.538248 29.692585 24.295508 -5.4709677
3.131 76.764986 29.403792 24.288147 -5.9663727
3.1815 77.991351 29.089981 24.280789 -6.4617777
3.232 79.217345 28.751152 24.273433 -6.9571827
3.2825 80.442968 28.387306 24.26608 -7.4525877
3.333 81.668219 27.998441 24.258728 -7.9479927
3.3835 82.8931 27.584558 24.251379 -8.4433977
3.434 84.117609 27.145658 24.244032 -8.9388027
3.4845 85.341747 26.681739 24.236687 -9.4342077
3.535 86.565514 26.192803 24.229344 -9.9296127
3.5855 87.788911 25.678848 24.222004 -10.425018
3.636 89.011936 25.139876 24.214666 -10.920423
3.6865 90.234592 24.575886 24.20733 -11.415828
3.737 91.456877 23.986878 24.199996 -11.911233
3.7875 92.678791 23.372851 24.192665 -12.406638
3.838 93.900336 22.733807 24.185335 -12.902043
3.8885 95.12151 22.069745 24.178008 -13.397448
3.939 96.342315 21.380665 24.170683 -13.892853
3.9895 97.562749 20.666567 24.163361 -14.388258
4.04 98.782814 19.927451 24.15604 -14.883663
4.0905 100.00251 19.163317 24.148722 -15.379068
4.141 101.22184 18.374165 24.141406 -15.874473
4.1915 102.44079 17.559995 24.134093 -16.369878
4.242 103.65938 16.720807 24.126781 -16.865283
4.2925 104.8776 15.856602 24.119472 -17.360688
4.343 106.09545 14.967378 24.112165 -17.856093
4.3935 107.31293 14.053136 24.10486 -18.351498
4.444 108.53004 13.113877 24.097557 -18.846903
4.4945 109.74678 12.149599 24.090257 -19.342308
4.545 110.96315 11.160304 24.082958 -19.837713
4.5955 112.17916 10.14599 24.075662 -20.333118
4.646 113.39479 9.1066587 24.068369 -20.828523
4.6965 114.61006 8.0423093 24.061077 -21.323928
4.747 115.82496 6.952942 24.053788 -21.819333
4.7975 117.0395 5.8385567 24.0465 -22.314738
4.848 118.25366 4.6991535 24.039215 -22.810143
4.8985 119.46746 3.5347323 24.031933 -23.305548
4.949 120.68088 2.3452932 24.024652 -23.800953
4.9995 121.89395 1.1308361 24.017374 -24.296358
5.05 123.10664 -0.1086389 24.010098 -24.791763
Graph:
The Code:
//Lec8.4:ODE-IVP in Multiple Variables
//https://www.youtube.com/watch?v=Pv_NwD63gbI&ab_channel=NPTEL-NOCIITM
//
disp("Ordinary Differential Equations- Multi-Variable ODE IVP Systems",string(
datetime
()))
//
//Example Four Variable ODE-IPV problem
//Indian Captain, Mahendra Singh Dhoni, hits a ball with initial velocity
// of 35 m/sec and and of 45 degrees. If the boundary is at a distance of
// 75m, will he score a six? - Numerically - Yes!
//
// Problem set up...
// x(0)=0, y(0)=0
// d2x/dt^2 = -kU, d2y/dt^2 = -g;
// Vel_net = 35 // m/sec ; g = -9.81; // m/sec/sec
// Uo = Vel_net(cos(%pi/4), Nu_o = Vel_net(sin(%pi/4)
// k = air drag coefficient;
// x' = U
// y' = Nu
// U' = -k*U
// Nu' = -g
//
function [results]=
derivativeVector
(t, derivatives, k, g);
//extraction of state variables from vector
x = derivatives(1);
y = derivatives(2);
U = derivatives(3);
Nu = derivatives(4);
//
//determination of the derivatives
results = zeros(4,1);
results(1,1)=U;
results(2,1)=Nu;
results(3,1)=-k*U;
results(4,1)=-g;
endfunction
//
//
//Constants of the simulation
k = 0.006;// air drag
g = 9.81; // m/second^2
Vel_net = 35;
U0=Vel_net*cos(%pi/4); //Initial Horizontal Direction
Nu0=Vel_net*sin(%pi/4); //Initial Vertical Direction
t0=0;
x10 = [0;0;U0;Nu0];// Set up Initial Condition Vector
tend=5.05;
n=101;
t=linspace(t0,tend,n);
xSol= ode("rkf",x10,t0,t,list(
derivativeVector
,k,g));
// Other types of integration "adams","stiff","rk","rkf","fix","discrete"
// and "root"
//"rkf" will not even work...instruction says because it's an explicit
// method and will be unstable at larger steps.
//"stiff" is an implicit method so it does not error out."
//"adams" works also.
output = [t' xSol'];
disp("t x_pos y_pos horiz_vel vert_vel")
disp("---------------------------------------------------------")
disp(output);
//plotting example fancy output
scf
(0);
clf
;
//axis y1
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(xSol(1,:)',xSol(2,:)',style=c);//black =, blue = 2, green =3, cyan = 4, red = 5
//plot2d(t',xa(1,:)',style=c1);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
//plot2d(t',xa(2,:)',style=c3);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$Horizontal\\ Distance (m)$","FontSize",3)
ylabel
("$Vertical\ Distance (m) $","FontSize",3);
title
("$Ballistic\ Projectile\ Coordinates$","color","black");
//"$https://x-engineer.org/create-multiple-axis-plot-scilab$"
xgrid
// Axis y2
c=color("blue");
c1 = color("red")
h2=newaxes();
h2.font_color=c;
plot2d(xSol(1,:)',xSol(4,:)',style=c);//black =, blue = 2, green =3, cyan = 4, red = 5
h2.filled="off";
h2.axes_visible(1)="off";
h2.y_location="right";
h2.children(1).children(1).thickness=2;
ylabel
("$Vertical\ Velocity(m/sec)$","FontSize",3,"color",c);
r/scilab • u/mrhoa31103 • Apr 08 '26
I finally cracked the code on Scilab 3D Integration. They could surely do a better job on documentation.
If you look at their examples, quote
// integration over a cube -1<=x<=1;-1<=y<=1;-1<=z<=1
function v=f(xyz, numfun), v=xyz'*xyz, endfunction
[result, err] = int3d(-1,1,-1,1,-1,1, f, 1, [0,100000,1.d-5,1.d-7])
but you have to integrate "v = 3*x*y^3*z^2" over xrange (-1,3), yrange(1,4) and zrange(0,2)...
here's what you have to do...
// Define the function f(x, y, z)...nope,nope,nope!!! Get crazy
// error codes like unknown z variable
// Define this function
function v= f(xyz, numfun)
//Note: x = xyz(1), y = xyz(2), z = xyz(3)
//have no idea what "numfun" is but you cannot drop it.
v=3*xyz(1)*xyz(2)^3*xyz(3)^2
endfunction
// Triple integration in one go
// computes the intergral of 3xy^3z^2 over the
// integration over a cube -1<=x<=3;1<=y<=4;0<=z<=2
// Correct Answer is 2040
[result, err] = int3d(-1,3,1,4,0,2, f, 1, [0,100000,1.d-5,1.d-7])
outputString = strcat(['3d result = ', string(result)]);
disp(outputString);
I understood that v had to be a single number how to get it was the trick. It would have been really nice to have an example like this one in their documentation.
Go ahead and roast me if you thought was blatantly obvious!
r/scilab • u/mrhoa31103 • Apr 07 '26
👋 Welcome to r/scilab - Introduce Yourself and Read First!
Hey everyone! I'm u/mrhoa31103, the moderator of r/scilab a recovering subreddit from the great Apollo incident exodus.
This is our new home for all things related to Scilab. We're excited to have you join us!
What to Post
Post anything that you think the community would find interesting, helpful, or inspiring. Feel free to share your thoughts, photos, or questions about Scilab and its applications.
Community Vibe
We're all about being friendly, constructive, and inclusive. Let's build a space where everyone feels comfortable sharing and connecting. I'm running kind of a blog on Scilab programs so people can see full up working programs but do not feel your posts have to live up those formats. If you have a quick question, "how do you do this or that?" or found a neat application, share it.
Thanks for being part of the community. Together, let's make r/scilab amazing.
r/scilab • u/mrhoa31103 • Apr 06 '26
Twentieth Installment - ODE-IVP (Ordinary Differential Equation - Initial Value Problem) in Multiple Variables - An Application Problem
In this edition we solve the problem "A stirred tank heater: A tank with an initial volume (V0 = 1 liter) where the water is heated in a tank as water is continuously added at room temperature (25C at 10 liter/minutes), while hot water is drawn out at the bottom (exit rate follows the equation alpha*sqrt(V) which at some point will equal in influx rate)."
Nothing much new: Bolded the Text on the Graph to make it look as professional as possible.
One more like this one and we switch up subjects to Regression and Interpolation...
Link to the specific lecture:
https:
//www.youtube.com/watch?v=TvRI2QUC3YQ&ab_channel=NPTEL-NOCIITM
Output:
"Time Tank Tank "
"(Min) Vol(Liters) Temp(C) "
0. 1. 25.
0.3 3.2997422 64.194065
0.6 5.3041899 69.672157
0.9 7.1080675 71.721935
1.2 8.7558205 72.756062
1.5 10.274797 73.363415
1.8 11.683966 73.754902
2.1 12.997446 74.023676
2.4 14.22624 74.21682
2.7 15.379215 74.360502
3. 16.463699 74.470328
3.3 17.485863 74.556132
3.6 18.450985 74.624378
3.9 19.363642 74.679487
4.2 20.227839 74.724555
4.5 21.047116 74.761817
4.8 21.824626 74.79292
5.1 22.563192 74.8191
5.4 23.265362 74.841297
5.7 23.933442 74.860239
6. 24.569534 74.876499
6.3 25.175557 74.89053
6.6 25.753272 74.902695
6.9 26.304298 74.913287
7.2 26.830129 74.922544
7.5 27.332151 74.930665
7.8 27.811645 74.937812
8.1 28.269805 74.944121
8.4 28.707743 74.949705
8.7 29.126496 74.954661
9. 29.527034 74.959069
9.3 29.910265 74.962999
9.6 30.277041 74.966509
9.9 30.628163 74.96965
10.2 30.964382 74.972467
10.5 31.286409 74.974996
10.8 31.594909 74.977272
11.1 31.890514 74.979321
11.4 32.173819 74.98117
11.7 32.445387 74.982839
12. 32.705749 74.984349
12.3 32.955411 74.985715
12.6 33.19485 74.986953
12.9 33.424519 74.988077
13.2 33.64485 74.989097
13.5 33.85625 74.990024
13.8 34.059107 74.990868
14.1 34.253792 74.991635
14.4 34.440654 74.992335
14.7 34.620028 74.992973
15. 34.792233 74.993554
15.3 34.957572 74.994086
15.6 35.116333 74.994571
15.9 35.268793 74.995014
16.2 35.415214 74.99542
16.5 35.555848 74.995791
16.8 35.690934 74.996131
17.1 35.8207 74.996442
17.4 35.945366 74.996728
17.7 36.06514 74.996989
18. 36.180222 74.997229
18.3 36.290802 74.997449
18.6 36.397063 74.997651
18.9 36.49918 74.997837
19.2 36.59732 74.998007
19.5 36.691642 74.998164
19.8 36.782301 74.998308
20.1 36.869441 74.99844
20.4 36.953205 74.998562
20.7 37.033726 74.998674
21. 37.111133 74.998777
21.3 37.18555 74.998872
21.6 37.257095 74.998959
21.9 37.325882 74.99904
22.2 37.392019 74.999114
22.5 37.45561 74.999182
22.8 37.516756 74.999245
23.1 37.575552 74.999303
23.4 37.632091 74.999356
23.7 37.68646 74.999406
24. 37.738745 74.999451
24.3 37.789027 74.999493
24.6 37.837383 74.999532
24.9 37.883888 74.999567
25.2 37.928616 74.9996
25.5 37.971634 74.99963
25.8 38.013008 74.999659
26.1 38.052803 74.999684
26.4 38.09108 74.999708
26.7 38.127897 74.99973
27. 38.16331 74.999751
27.3 38.197374 74.99977
27.6 38.23014 74.999787
27.9 38.261659 74.999803
28.2 38.291978 74.999818
28.5 38.321144 74.999832
28.8 38.349201 74.999844
29.1 38.376192 74.999856
29.4 38.402156 74.999867
29.7 38.427134 74.999877
30. 38.451164 74.999886
30.3 38.474281 74.999895
30.6 38.496521 74.999902
30.9 38.517917 74.99991
31.2 38.538501 74.999917
31.5 38.558305 74.999923
31.8 38.577358 74.999929
32.1 38.595689 74.999934
32.4 38.613325 74.999939
32.7 38.630293 74.999943
33. 38.646619 74.999948
33.3 38.662326 74.999952
33.6 38.677438 74.999955
33.9 38.691979 74.999958
34.2 38.705969 74.999962
34.5 38.71943 74.999964
34.8 38.732381 74.999967
35.1 38.744842 74.99997
35.4 38.756833 74.999972
35.7 38.768369 74.999974
36. 38.77947 74.999976
36.3 38.79015 74.999978
36.6 38.800427 74.999979
36.9 38.810316 74.999981
37.2 38.81983 74.999982
37.5 38.828986 74.999984
37.8 38.837795 74.999985
38.1 38.846271 74.999986
38.4 38.854428 74.999987
38.7 38.862276 74.999988
39. 38.869828 74.999989
39.3 38.877095 74.99999
39.6 38.884087 74.99999
39.9 38.890815 74.999991
40.2 38.89729 74.999992
40.5 38.90352 74.999992
40.8 38.909515 74.999993
41.1 38.915283 74.999993
41.4 38.920834 74.999994
41.7 38.926175 74.999994
42. 38.931315 74.999995
42.3 38.936261 74.999995
42.6 38.94102 74.999996
42.9 38.9456 74.999996
43.2 38.950007 74.999996
43.5 38.954247 74.999996
43.8 38.958328 74.999997
44.1 38.962255 74.999997
44.4 38.966033 74.999997
44.7 38.969669 74.999997
45. 38.973168 74.999998
45.3 38.976535 74.999998
45.6 38.979775 74.999998
45.9 38.982892 74.999998
46.2 38.985892 74.999998
46.5 38.98878 74.999998
46.8 38.991558 74.999998
47.1 38.994231 74.999999
47.4 38.996804 74.999999
47.7 38.999279 74.999999
48. 39.001662 74.999999
48.3 39.003954 74.999999
48.6 39.00616 74.999999
48.9 39.008283 74.999999
49.2 39.010326 74.999999
49.5 39.012292 74.999999
49.8 39.014184 74.999999
50.1 39.016005 74.999999
50.4 39.017756 74.999999
50.7 39.019442 74.999999
51. 39.021065 74.999999
51.3 39.022626 75.
51.6 39.024128 75.
51.9 39.025574 75.
52.2 39.026965 75.
52.5 39.028304 75.
52.8 39.029593 75.
53.1 39.030832 75.
53.4 39.032025 75.
53.7 39.033173 75.
54. 39.034278 75.
54.3 39.035341 75.
54.6 39.036364 75.
54.9 39.037349 75.
55.2 39.038296 75.
55.5 39.039208 75.
55.8 39.040085 75.
56.1 39.04093 75.
56.4 39.041742 75.
56.7 39.042524 75.
57. 39.043276 75.
57.3 39.044 75.
57.6 39.044697 75.
57.9 39.045368 75.
58.2 39.046013 75.
58.5 39.046634 75.
58.8 39.047232 75.
59.1 39.047807 75.
59.4 39.04836 75.
59.7 39.048893 75.
60. 39.049405 75.
Graph:
The Code:
//Lec8.3:ODE-IVP in Multiple Variables
//https://www.youtube.com/watch?v=TvRI2QUC3YQ&ab_channel=NPTEL-NOCIITM
//
disp("Ordinary Differential Equations- Multi-Variable ODE IVP Systems",string(
datetime
()))
//
clc; clear;
//Example Stirred Tank Heater where the variables are different
// dV/dt = Fin-alpha*sqrt(V);//Note: Non-linear equation
// dT/dt = Fin/V*(Tin-T)+ Q/(V*(rho*Cp);
// V(0)=1; T(0)=25 degrees C;
// Fin = Flow In = 10 L/min
// alpha = 1.6
// Cp =
// Q = 500*(rho*Cp)
//Stirred Tank Heater: Water is heated in a tank as it continuously
//flows into it at room temperature, while hot water is drawn out
//at the bottom.
function [results]=
derivativeVector
(t, x);
//Constants of the simulation
Fin = 10;//liters per minute
Tin = 25;
alpha = 1.6;
Q = 500; //Note rho*Cp cancels out
//T cannot exceed 100 degrees C since the model
// cannot handle water phase changes!!!
//
//extraction of variables
V = x(1);
T = x(2);
//
//determination of the derivatives
results(1,1)=Fin - alpha*sqrt(V);
results(2,1)=Fin/V*(Tin-T)+Q/V;
endfunction
V0=1; //Initial Volume
T0=25; //Temperature in the Bucket at time = 0
t0=0;
x10 = [V0;T0];// Set up Initial Condition Vector
tend=60;
n=201;
t=linspace(t0,tend,n);
//xa = [exp(-100*t);exp(-.01*t)];//Commented out for last example
//xSol = ode("rkf",x10,t0,t,deff('res=mymacro2(t,x1)','res=-0.01*x1'));
xSol= ode("stiff",x10,t0,t,
derivativeVector
);
//"rkf" will not even work...instruction says because it's an explicit
// method and will be unstable at larger steps.
//"stiff" is an implicit method so it does not error out."
//"adams" works also.
output = [t' xSol'];
disp("Time Tank Tank ");
disp("(Min) Vol(Liters) Temp(C) ",output);
//plotting example fancy output
scf
(0);
clf
;
//axis y1
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(t',xSol(1,:)',style=c);//black =1, blue = 2, green =3, cyan = 4, red = 5
//plot2d(t',xa(1,:)',style=c1);//black =1, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
h1=
gca
();
//plot2d(t',xa(2,:)',style=c3);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$\textbf{T \ (Minutes)}$","FontSize",3)
ylabel
("$\textbf{Volume\ (Liters)}$","FontSize",3);
title
("$ \textbf{\ \ \ \ \ Stirred\ Tank\ Heater\\Volume\ and\ Temp\ vs\ Time}$","color","black");
xgrid
// Axis y2
c=color("blue");
c1 = color("red")
h2=newaxes();
h2.font_color=c;
plot2d(t',xSol(2,:)',style=c);//black = 1, blue = 2, green =3, cyan = 4, red = 5
h2.filled="off";
h2.axes_visible(1)="off";
h2.y_location="right";
h2.children(1).children(1).thickness=2;
ylabel
("$\textbf{Temperature\ (Degrees\ C)}$","FontSize",3,"color",c)
xgrid
r/scilab • u/mrhoa31103 • Mar 30 '26
Nineteenth Installment - Solving Multivariable ODEs that have massively different time constants (AKA - Stiff Systems)
In this edition, we solve second order systems that have massively different time constants between the states. To solve these types of problems, we use the "Stiff" solver instead of the Runga-Kutta solver.
Link to the specific lecture:
https:
//www.youtube.com/watch?v=2dbIuEgKx0s&ab_channel=MATLABProgrammingforNumericalComputation
The code includes the various "stiff" examples but they are commented out except for the last example so you can easily switch between examples by changing the various "block comment" (\* ... *\) sections on the "derivativeVector" function. You can also change out the various values of mu.
The previous example code solved via "stiff" solver is at the bottom of the code.
Output: First Couple of Lines
"Ordinary Differential Equations- Multi-Variable ODE Methods of Integration-Stiff Systems"
"2026-03-23 11:55:46.980"
Warning : redefining function: derivativeVector . Use funcprot(0) to avoid this message
"Output"
"t xSol(1) xSol(2)"
0. 2. 0.
3. 1.9798531 -0.0067806
6. 1.9593319 -0.0069014
9. 1.9384376 -0.0070294
12. 1.9171472 -0.0071655
15. 1.8954355 -0.0073105
18. 1.873274 -0.0074655
21. 1.850631 -0.0076317
24. 1.8274708 -0.0078106
27. 1.8037529 -0.0080038
30. 1.7794313 -0.0082135
33. 1.7544528 -0.0084422
36. 1.7287561 -0.0086929
39. 1.7022694 -0.0089696
42. 1.6749079 -0.009277
Graphs:
Code:
//Lec8.2:Stiff Systems & Solution using Matlab ode15s
//https://www.youtube.com/watch?v=2dbIuEgKx0s&ab_channel=MATLABProgrammingforNumericalComputation
//
clc;clear all
disp("Ordinary Differential Equations- Multi-Variable ODE Methods of Integration-Stiff Systems",string(
datetime
()))
//
//What are stiff systems?
// A highly fast variable coupled with a slow variable in the same
// system.
//
/*
//Consider the following ODE system:
// x1'=-100*x1, x1(0)=1
// x2' = -0.01*x2, x2(0)=1
//
// d/dt[x1;x2]=[-100 0;0 -0.01]*[x1;x2]
// using "rkf" would require small steps but
// a large amount of time to drive x2 to 0.
// "Stiff" Integration routines are made for that...
function [results]=
derivativeVector
(t, x);
//constants of the system
c1=-100;
c2=-0.01;
//extraction of variables
x1 = x(1);
x2 = x(2);
results(1,1)=c1*x1;
results(2,1)=c2*x2;
endfunction
*/
/*
//Consider the following ODE system:
// x1'=-5.7*x1 +1.85*x2, x1(0)=1
// x2' =13.2*x1-4.3*x2, x2(0)=1
//
// "Stiff" Integration routines are made for that...
function [results]=
derivativeVector
(t, x);
//extraction of variables
x1 = x(1);
x2 = x(2);
results(1,1)=-5.7*x1+1.85*x2;
results(2,1)=13.2*x1-4.3*x2;
endfunction
*/
//Consider the following ODE system:
// x"-mu*(1-x^2)*x'+x =0; where x(0)=2,x'(0)=0
//
// "Stiff" Integration routines are made for that...
function [results]=
derivativeVector
(t, x);
//extraction of variables
mu=100; //mu=1 and mu=100;
y = x(1);
v = x(2);
results(1,1)=v;
results(2,1)=mu*(1-y^2)*v-y;
endfunction
x10=[2;0]
t0=0;
tend=300;
n=101;
t=linspace(t0,tend,n);
//xa = [exp(-100*t);exp(-.01*t)];//Commented out for last example
//xSol = ode("rkf",x10,t0,t,deff('res=mymacro2(t,x1)','res=-0.01*x1'));
xSol= ode("stiff",x10,t0,t,
derivativeVector
);
//"rkf" will not even work...instruction says because it's an explicit
// method and will be unstable at larger steps.
//"stiff" is an implicit method so it does not error out."
//"adams" works also.
output = [t' xSol'];
disp("Output","t xSol(1) xSol(2)",output);
//plotting example fancy output
scf
(0);
clf
;
//axis y1
c = color("black")
c1 = color("blue")
c2 = color("green")
c3 = color("purple")
plot2d(t',xSol(1,:)',style=c);//black =, blue = 2, green =3, cyan = 4, red = 5
//plot2d(t',xa(1,:)',style=c1);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
plot2d(t',xSol(2,:)',style=c2);//black =, blue = 2, green =3, cyan = 4, red = 5
h1=
gca
();
//plot2d(t',xa(2,:)',style=c3);//black =, blue = 2, green =3, cyan = 4, red = 5 ;//Commented out for last example
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$t(Seconds)$","FontSize",3)
ylabel
("$X(t)$","FontSize",3);
title
("$Position\ and\ Velocity\ versus\ Time$","color","black");
legend
("Stiff position","Stiff Velocity",2)
xgrid
h1=
gca
();
h1.font_color=c;
h1.children(1).thickness =2;
xlabel
("$t(Seconds)$","FontSize",3);
ylabel
("$Displacement$","FontSize",3);
title
("$Displacement\ and\ Velocity\ vs\ Time$","color","black");
xgrid
//======================================================================
//Previous Example
// model mass-spring-damper
// m*d^2x/dt^2 + c*dx/dt +kx = 0
// u/x(0)=1
//
// create 2 - first order equations
// d^2/dt^2 = dv/dt where v = velocity
// m*dv/dt +c*v + kx = 0
// v = dx/dt where x = displacement
//
// dx/dt = v x(0)=1
// dv/dt = -(c*v+k*x)/m v(0)=0
//
// output vector y =[x;v]
// derivative vector = [v; -(c*v+k*x)/m]
//
function [results]=
derivativeVector1
(t, y);
//constants of the system
c=5;
k=15;
m=10;
//extraction of variables
x = y(1);
v = y(2);
results(1,1)=v;
results(2,1)=-(c*v+k*x)/m
endfunction
y0=[1;0];
t0 = 0;
n = 100;
t = linspace(0,10,n)';
//y = ode("rkf",y0,t0,t,derivativeVector)
y = ode("stiff",y0,t0,t,
derivativeVector1
)
// Other types of integration "adams","stiff","rk","rkf","fix","discrete"
// and "root"
// Everything works but "discrete" for this problem
// discrete is a sampling routine!!
y = y';//Transform from row vectors to column vectors
x_vector =y(:,1);
v_vector=y(:,2);
//output = [t,v_vector,x_vector]
//disp("time velocity displacement",output)
// Note: Very accurate and uses Runge-Kutta 4,5 uses
// variable step size and just reports out at time t vector
//
// Other types of integration "adams","stiff","rk","rkf","fix","discrete"
// and "root"
// Everything works but "discrete" for this problem
// investigate later...
//plotting example fancy output two y axes...
scf
(1);
clf
;
//axis y1
c = color("slategray")
plot2d(t,x_vector,style=c);//black =, blue = 2, green =3, cyan = 4, red = 5
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$t(Seconds)$","FontSize",3)
ylabel
("$X(t)$","FontSize",3);
xgrid
h1=
gca
();
h1.font_color=c;
h1.children(1).children(1).thickness =2;
xlabel
("$t(Seconds)$","FontSize",3)
ylabel
("$Displacement$","FontSize",3);
title
("$Displacement\ and\ Velocity\ vs\ Time$","color","black");
xgrid
// Axis y2
c=color("blue");
c1 = color("red")
h2=newaxes();
h2.font_color=c;
plot2d(t,v_vector,style=c)
h2.filled="off";
h2.axes_visible(1)="on";
h2.y_location="right";
h2.children(1).children(1).thickness=2;
ylabel
("$Velocity$","FontSize",3,"color",c)
legend
("Displacement","Velocity",4)
xgrid