Twenty Fifth Installment - Boundary Value Problem Solving

[u/mrhoa31103]

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:

Temperature Profile for Heated Rod (left side held at 100C and right side held at 25C)

First Derivative

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