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

2 views (last 30 days)
SARSKOLIN FOSSO on 26 Sep 2022
Commented: SARSKOLIN FOSSO on 29 Sep 2022
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
TT=T(2:n,k-1);
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;
end
end
figure(1)
plot(t,u)
grid figure(2)
plot(t,Tq)
grid figure(3)
surf(x,t,T) .
SARSKOLIN FOSSO on 28 Sep 2022
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

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?
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)