Coupled first order differential equations

[Rafique Mir]

How i can solve a system of 6 first order differential equations by using numerical techniques like Euler method, RK-4th order method , ODE -45 etc.How i can solve a system of 6 first order differential equations by using numerical techniques like Euler method, RK-4th order method , ODE -45 etc.

 

[mathmike]

can you give an example

 

[Rafique Mir]

Coupled ODEs

Dear here are the system of six differential equtions i want to solve them by Euler method , Predictor corrector method and for RK-4th order method. I have made a program in Matlab please check wether it is right or not.
function [GG ] = myfun_rafiq(t,x)

% ----- State Variables Selection
% x1 = nw, x2 = Nw, x3 = nc, x4 = Nc, x5 = np, x6 = Np

% ----- Parameters

sigmaEjQj = 793.7;
EiQi =700;
EpQp = 13.7;
EcQc = 80;
EfQf = 0;
Vc = 9.08e6;
Vw = 1.37e7;
Vp = 1.37e6;
lk = 0.0;
Kc = 40.0;
Kp = 6.90;
f_t = 1;
sigma = 13.4e-24;
fn = 1;
fs = 0.5;
phi0 = 9.2e13;
phiE = 0.0026*phi0;
lambda = 7.4612e-005; % PER SECOND

N0 = 6.023e23;
A = 56;
S = 1.01e8;

% ---- parameters + definitio of C(t)

a = 50*3600; % in seconds
deltaCDeltaT = 20e-12;
% deltaCDeltaT=deltaCDeltaT/(3600^2);
C0 = 2.4e-13;
Cs = 25e-6;
b = 400*3600; % in seconds
t0 =50*3600; % in seconds %REPLCED 500 with 50
% t = 900*3600;
if t < a
C_t = 0;
elseif (t > a & t < b)
C_t = deltaCDeltaT*(t-t0);
else
C_t = Cs;
end
% t
% C_t

% ---- Parameters + definition of g(t)
%
tin = 500*3600;
tmax = 600*3600;
w0 = 18.3e6;
w2 = 0.3*w0;
alpha = 0.004;
% g_t = 1;

if t <= tin
g_t = 1;
elseif (t > tin & t <= tmax)
g_t = 1 - alpha*(tin - t);
else %REPLACED 'b' with tmax
g_t = w2/w0;
end
% if(t<tin)
% g_t = 1;
% elseif((t<=tmax))
% g_t = 1 - alpha*(tin - t);
% % g_t=0.00001;
% else %REPLACED 'b' with tmax
% g_t = w2/w0;
% % g_t=0.00001;
% end

% % pause

Sw = (C_t* S*N0*fn*fs)/(Vw*A);

dxdt = zeros(6,1);

dxdt(1) = sigma*phiE*x(2) - ( (sigmaEjQj*g_t/Vw) + lk*g_t/Vw + lambda)*x(1) + (Kp*g_t/Vw)*x(5) + (Kc*g_t/Vw)*x(3);

dxdt(2) = -( (sigmaEjQj*g_t/Vw) + lk*g_t/Vw + sigma*phiE)*x(2) + (Kp*g_t/Vw)*x(6) + (Kc*g_t/Vw)*x(4) + Sw;

dxdt(3) = sigma*phi0*x(4) + (EcQc*g_t/Vc)*x(1) - ( Kc*g_t/Vc + lambda)*x(3);

dxdt(4) = (EcQc*g_t/Vc)*x(2) - ( Kc*g_t/Vc + sigma*phi0)*x(4);

dxdt(5) = (EpQp*g_t/Vp)*x(1) - ( Kp*g_t/Vp + lambda)*x(5);

dxdt(6) = (EpQp*g_t/Vp)*x(2) - ( Kp*g_t/Vp)*x(6);

GG=[dxdt(1); dxdt(2); dxdt(3); dxdt(4); dxdt(5); dxdt(6)];

% ---- PIEAS

% Clearing work space

clc;
clear all;
% close all;
tic
% Define time for simulation
%x0 = [0.1, 0.2, 0.3, 0.4, 0.1, 0.1]*10^-7;
% x0=[0 0.5 0 0.5 0 0.5];
% T = 1:500*3600; % 10 seconds
% T=T*3600;

% Define initional conditions

%R=myfun_rahat(T,x0);

% Run simulation
% options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);

[tt,xx] = ode45(@myfun_rafiq,[0 1000*3600],[0 0.5 0 0.5 0 50]);

toc

% Plot Results
grid on;
figure(1),plot( tt, xx(:,1),'r.-',tt,xx(:,3),'m:',tt,xx(:,5),'b.-');
xlabel('t(sec)');
title('Specific Activity');
legend('nw','nc','np');
figure(2)
grid on;
plot(tt, xx(:,1),'r-');
%axis([0 1.81e5 0 3e-4]);

 

[saltydog]

Dear here are the system of six differential equtions i want to solve them by Euler method , Predictor corrector method and for RK-4th order method. I have made a program in Matlab please check wether it is right or not.

dxdt(1) = sigma*phiE*x(2) - ( (sigmaEjQj*g_t/Vw) + lk*g_t/Vw + lambda)*x(1) + (Kp*g_t/Vw)*x(5) + (Kc*g_t/Vw)*x(3);

Oh my God dude! Now, I really admire you for taking the time to input all that code but it's just too hard to read. I picked out the first one above and I can't understand it. I realize you're probably not familiar with the math-formatting code we use in here called "LaTex" but that would help with understanding your code. If you want to try using it, you can go to the General Physics forum and read the "Introducing LaTex" thread at the top. Also, I use Mathematica so would not be able to help you with Matlab.

 

[HallsofIvy]

If you are already familiar with numerical methods for a single equation in one variable, just repeat them: For 6 coupled equations in 6 unknowns, set up 6 Runge-Kutta iterations, doing all 6 at each step then using the results from all 6 for the values of the 6 variables in the next iteration.

 

[peterbone]

I had exactly the same problem when trying to simulate a double inverted pendulum a few months ago. The system is 3 coupled second order differential equations which I reduced to 6 first order differential equations. Here's the program (written in Delphi).
http://atlas.walagata.com/w/peterbone/Balance.zip
The code for the double inverted pendulum is in the DoublePendulum.pas file which you can open in notepad. You can then see how the equations are solved. You should be able to understand it even if you don't know pascal.