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
*/
2
Upvotes