please give me your suggestions. thanks

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

Show older comments

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)

.

##### 3 Comments

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 ?

### Answers (1)

Bill Greene
on 28 Sep 2022

##### 3 Comments

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

### See Also

### Community Treasure Hunt

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

Start Hunting!