Hello everyone I hope you are all doing well. Please I have a concern about the coupling of an ODE and an EDP.

1 view (last 30 days)
I would like to couple a second order differential equation (ODE) with the heat equation (PDE). For this, I used the RK4 method for the ODE and the Crank-Nicholson methode for the EDP. I don't know why my code does not give me a good results. Please give me a hand. I have attached two files to this request: one which presents my code and the other which shows the equation system
close all
clear all
% Matlab code : The dynamic of slider using the heat equation
% I would like to use the Cranck Nicholson methode to solve the heat equation and
% incorporate it simultaneously into the second order of ODE
% Parameters to define the dynamic of slider.
Kr = 10; m = 10; F0 = 1; nu0 = 10; Vp = 10.^(-11); gama = 0.8; beta = 0.2; uc = 0.5;
% Parameters to define the heat equation.
rho = 2700; cs = 1000; Kc = 2.7;
% Initial temperature of heat
Ti = 150;
% boundary condition of heat equation.
Tc = 1050; Th = 1150;
% final lenght and time.
L = 2; % thickness of conduction chamber.
tf =80; % Final time.
% Parameters needed to solve the heat equation.
maxk = 100; % Number of time steps
dt = tf/maxk; % lenght on time steps
n = 100; % Number of space steps
dx = L/n; % lenght on space steps
% Definition of time
t = 0:dt:maxk.*dt;
% parameters used to make the reduction of another parameters.
H = L./2; % the half thickness of conduction chamber
Kd = Kc./(rho.*cs); % thermal diffusivity
alpha = (Kd.*dt)./((dx).^(2)); % reduction parameter
alpha1 = (alpha*dx*dx)./(Kc*H); % reduction parameter
% % Previsional matrice to incorporate data.
u = zeros(1,maxk+1);
v = zeros(1,maxk+1);
T = zeros(n+1,maxk+1);
% Defining the Matrices M_left :(ML)
aal(1:n-2) = -alpha;
bbl(1:n-1) = 2 + 2.*alpha;
ccl(1:n-2) = -alpha;
ML = diag(bbl,0)+diag(aal,-1)+diag(ccl,1);
% Defining the Matrices M_rigth :(ML)
aar(1:n-2) = alpha;
bbr(1:n-1) = 2 - 2.*alpha;
ccr(1:n-2) = alpha;
MR = diag(bbr,0)+diag(aar,-1)+diag(ccr,1);
% initial slip, velocity and enery.
u(1) = 0; % slip.
v(1) = 0; % velocity.
Tq(1) = Th; % temperature
% Definition of the three functions using by Rk4.
pascak1 = @(t,u,v) v;
pascak2 = @(t,u,v) (-1./m).*((Kr-gama).*u + nu0.*exp(beta*((Th-Tq)./(Th-Tc))).*v + F0 - Kr.*Vp.*t);
for i = 1:n+1
% Initial temperature of the fault.
x(i) =(i-1)*dx;
T(i,1) = Ti;
% Boundary condition
for k=2:maxk + 1 % Time Loop
T(1,k) = Th;
T(n+1,k) = Tc;
% time(k) = (k-1)*dt;
L1(k) = dt.*pascak1(t(k-1),u(k-1),v(k-1));
K1(k) = dt.*pascak2(t(k-1),u(k-1),v(k-1));
L2(k) = dt.*pascak1(t(k-1) + dt/2, u(k-1) + L1(k)./2 ,v(k-1) + K1(k)./2);
K2(k) = dt.*pascak2(t(k-1) + dt/2, u(k-1) + L1(k)./2 ,v(k-1) + K1(k)./2);
L3(k) = dt.*pascak1(t(k-1) + dt/2, u(k-1) + L2(k)./2 ,v(k-1) + K2(k)./2);
K3(k) = dt.*pascak2(t(k-1) + dt/2, u(k-1) + L2(k)./2 ,v(k-1) + K2(k)./2);
L4(k) = dt.*pascak1(t(k-1) + dt, u(k-1) + L3(k) ,v(k-1) + K3(k));
K4(k) = dt.*pascak2(t(k-1) + dt, u(k-1) + L3(k) ,v(k-1) + K3(k));
% Implementation of Crank-Nicholson method
T(2:n,k)= inv(ML)*MR*TT + alpha1*F0*inv(ML)*v(k-1)*ones(99,1) - alpha1*gama*inv(ML)*u(k-1)*v(k-1)*ones(99,1);
% Because T is two dimenssions I extract Tq in one dimenssion at a specific space and
% incorporate it in the ODE equation.
Tq = T(50,1:maxk+1);
u(k) = u(k-1) + (1/6).*(L1(k) + 2*L2(k) + 2*L3(k) + L4(k));
v(k) = v(k-1) + (1/6).*(K1(k) + 2*K2(k) + 2*K3(k) + K4(k));
t(k) = t(k-1) + dt;
Torsten on 26 Sep 2022
Edited: Torsten on 26 Sep 2022
My suggestion would be that you don't use your own Runge-Kutta solver for the differential equations, but MATLAB's ODE15s. This will prevent at least one source of error.
And I don't know why you get one-dimensional arrays as solutions for u, v and T. Since your problem depends on the space and the time coordinate, they should be two-dimensional, shouldn't they ?
SARSKOLIN FOSSO on 28 Sep 2022
Thank you very much for your answer Professor.
I have one-dimensional tables for u and v because I extract the values ​​of T in the form of scalars (Tq as in my code) and I gradually put them back in the equation (1) and (2) and vice versa.
But it doesn't do it gradually. maybe it could be the execution sequence of the for loops.
And I don't know how to do.
I have attached a file which shows the behaviour of my system.
thank a lot for your help

Sign in to comment.

Answers (1)

Bill Greene
Bill Greene on 28 Sep 2022
I looked at your mathematical description of this problem.Since T is a function of z, V and u are also functions of z. So you really are solving three PDE. I see that in your code, you use T at the mid-point of the region in the second equation to convert it to an ODE. Are you sure this doesn't introduce significant errors?
Torsten on 28 Sep 2022
Since T appears in the equation for V, V becomes automatically two-dimensional (and so U).
Thus U = U(t,z) and V = V(t,z).
Otherwise I do not know which value for T you would insert in the equation for V (maybe a mean value of T over the spatial coordinate z ?).
SARSKOLIN FOSSO on 29 Sep 2022
Thank you professor for this question.
In my code I take in T(z,t) a scalar Tq(50,t) that I put back in the two other equations.
this is why I have u, v and Tq as scalars. I thus leave from two dimensions to one dimension.
my real problem is that I don't know how to properly sequence the execution.
thank you for your suggestion that I can also turn all these equations in two dimenssion.
now I am going to find how can I obtain U(z=0,t), V(z=0,t) and T(z,t)

Sign in to comment.


Find more on Programming in Help Center and File Exchange


Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!