# Solving Finite Difference Method Using ODE15s

19 views (last 30 days)
Mathushan Sabanayagam on 6 Aug 2019
Hi,
I have created a function in which I am solving a partial differential equation where temperature is dependent on time and radius (energy balance in spherical coordinates). I have discretised the spatial coordinates into n nodes. This yielded an ode with respect to time that I must complete over each node. There are also two other differential equations (not relevant to this however, its solutions are stored in n+1 and n+2 of the DyDt matrix just for information). I have shown the relevant parts of the code for brevity:
T = zeros(n,1); %initialise T as a matrix
DTDt = zeros(n,1); %initialise DTdt as a matrix
DyDt = zeros(n+2,1); %initialise a matrix containing DTdt from 1:n. (the n+1 and n+2 are the two other solutions)
T = y(1:n); %fill the T values into a y matrix from 1 to nth column
I then have a for loop for i=1:n-1 with my expression for DTDt. This allows me to solve all nodes up unti n-1. However, my problem is for the solution to T(n), I do not have a DTdt equation, but instead I have an algebraic equation:
T(n)=((2.*h.*Tg.*dr)-(ks.*T(n-2))+(4.*ks.*T(n-1)))/((3.*ks)+(2.*h.*dr));
It relies only only values from the other nodal positions and some other constants.
In a separate file, I utilise the functions, setting the initial condition as T0 at 298K.
How can I go about solving such a system? I understand it is a set of differential equations and 1 algebraic equation.
Thank you very much in advance. I am a beginner in Matlab and will appreciate any help.

Torsten on 8 Aug 2019
Edited: Torsten on 8 Aug 2019
%OTHER FILE THAT CALLS FUNCTION (Excluding all the constants)
T0 = ones(n,1)*298; %matrix of initial condition for T, rows = number of nodes, only 1 column since t=0
a0 = 1; %initial condition for a
b0 = 1; %initial condition for b
y0 = [T0;a0;b0]; %T0 followed by a0 and b0 in a single column matrix
M = eye(n+2);
M(n,n) = 0.0;
options = odeset('Mass',M)
%Call the solver
[T Y] = ode15s(@(t,y)convection(t,y,n,preexpomelt,activationmelt,preexpodecom,enthalpymelt,enthalpydecom,activationdecom,heatrate,alphas,alphal,dr,cps,R0,cpl,Tm,h,Tg,ks),t,y0,options)
And in "convection", set
dTdt(n) = T(n)-(((2.*h.*Tg.*dr)-(ks.*T(n-2))+(4.*ks.*T(n-1)))/((3.*ks)+(2.*h.*dr)));
This is the second method which should be easier to implement since the solver still has the same number of unknowns as before (n+2).
Mathushan Sabanayagam on 9 Aug 2019
Thank you for explaining, it makes sense. Thank you for all your help, really appreciate it!

Torsten on 7 Aug 2019
Edited: Torsten on 7 Aug 2019
You can either solve for T1,...,T_(n-1) (and the two extra variables) using ODE15S and calculate
T(n) internally each time the function to calculate the time derivatives is called or you can use the mass matrix option of ODE15S by setting M(n,n) = 0 instead of M(n,n) = 1. This will tell ODE15S that equation n is an algebraic equation instead of a differential equation.
Mathushan Sabanayagam on 7 Aug 2019
Edited: Mathushan Sabanayagam on 7 Aug 2019
Hi Torsten,
Thank you for your reply! Sorry if its an obvious question but how would I go about calculating the T(n) internally. Would the expression be in the function file or the file that calls the function? I would prefer the first method as I am already using ODE15s and so I would find it easier to implement. Thanks in advance for your help!
I've attached the code if that helps:
function DyDt=convection(t,y,n,preexpomelt,activationmelt,preexpodecom,enthalpymelt,enthalpydecom,activationdecom,heatrate,alphas,alphal,dr,cps,R0,cpl,Tm,h,Tg,ks)
%INITIALISING VARIABLES AND DERIVATIVES AS MATRICES
T = zeros(n,1);
a=0;
b=0;
DTDt = zeros(n,1);
DbDt=0;
DyDt = zeros(n+2,1);
DTDr = zeros(n-1,1);
D2TDr2 = zeros(n-1,1);
%STORING SOLUTIONS IN Y MATRIX
T = y(1:n); %fill the T values into a y matrix from 1 to nth element
a = y(n+1); %fill the a values into a y matrix at the n+1th element
b= y(n+2); %fill the a values into a y matrix at the n+2th element
for i=2:n-1
DTDr(i) = ( T(i+1) - T(i) )./ (dr) ;
D2TDr2(i) = ( 1./(dr^2) ) .* ( T(i+1) - (2.*T(i)) + T(i-1) );
end
%TEMPERATURE TIME DERIVATIVE
for i=2:n-1
DTDt(i) = ((((((a.^3)/(b.^3)).*alphas) + ((((b.^3)-(a.^3))/b.^3).*alphal))./(((R0.*b).^2))).* (D2TDr2(i))) + ((((((a.^3)./(b.^3)).*alphas) + ((((b.^3)-(a.^3))./b.^3).*alphal))./(((R0.*b).^2))).* (DTDr(i).*(2./(i.*dr))))...
+ (((1./((((a.^3)./(b.^3)).*cps) +((((b.^3)-(a.^3))./b.^3).*cpl)))).*((enthalpymelt.*preexpomelt.*exp(-activationmelt./(8.314.*T(i))))+(enthalpydecom.*preexpodecom.*exp(-activationdecom./(8.314.*T(i))))));
end
for i=1
DTDt(i)=DTDt(i+1);
end
T(n)=((2.*h.*Tg.*dr)-(ks.*T(n-2))+(4.*ks.*T(n-1)))/((3.*ks)+(2.*h.*dr));
if T(n) < Tm
else
end
DbDt(1)= ((-1.*b.*preexpodecom)./3).*exp(-activationdecom./((8.314.*T(n))));
end
%OTHER FILE THAT CALLS FUNCTION (Excluding all the constants)
T0 = ones(n-1,1)*298; %matrix of initial condition for T, rows = number of nodes, only 1 column since t=0
a0 = 1; %initial condition for a
b0 = 1; %initial condition for b
y0 = [T0;a0;b0]; %T0 followed by a0 and b0 in a single column matrix
%Call the solver
[T Y] = ode15s(@(t,y)convection(t,y,n,preexpomelt,activationmelt,preexpodecom,enthalpymelt,enthalpydecom,activationdecom,heatrate,alphas,alphal,dr,cps,R0,cpl,Tm,h,Tg,ks),t,y0)