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

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]

please give me your suggestions. thanks

[Torsten]

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]

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

Answer 1

[Bill Greene]

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]

Thank you very much for the time you are taking for this problem.

with each question (Torsten and Bill Greene) I think more.

U(t) and V(t) are one-dimensional (depending only on t) it is only T(z,t) which depends on two dimensions. Thank you for this remark I don't really know if the valuesT(50,t) which are scalars gradually enter the equation 1 and 2 and vice versa. It is an important question to ask

I have attached a file (picture) which shows the behaviour of my system (u(t) and v(t) one dimenssion; T(z,t) two dimenssions

[Torsten]

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]

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)