- #1

Rafique Mir

- 5

- 0

You are using an out of date browser. It may not display this or other websites correctly.

You should upgrade or use an alternative browser.

You should upgrade or use an alternative browser.

- Thread starter Rafique Mir
- Start date

- #1

Rafique Mir

- 5

- 0

- #2

mathmike

- 208

- 0

can you give an example

- #3

Rafique Mir

- 5

- 0

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 plz 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]);

- #4

saltydog

Science Advisor

Homework Helper

- 1,591

- 3

Rafique Mir said: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 plz 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.

- #5

HallsofIvy

Science Advisor

Homework Helper

- 43,021

- 970

- #6

peterbone

- 13

- 0

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.

Share:

- Replies
- 1

- Views
- 2K

- Replies
- 2

- Views
- 7K

- Replies
- 11

- Views
- 39K

- Last Post

- Replies
- 4

- Views
- 3K

- Last Post

- Replies
- 11

- Views
- 1K

- Last Post

- Replies
- 1

- Views
- 460

- Last Post

- Replies
- 2

- Views
- 2K

- Replies
- 1

- Views
- 3K

- Last Post

- Replies
- 4

- Views
- 1K

- Last Post

- Replies
- 2

- Views
- 2K