0

I am trying to solve 3 differentail equations(Lorenz equations) using ode solver: ode23s in Matlab. Here are the 3 lorenz equations:

dc/dt= alpha*I*(1-c) + c*(- k_f - k_d - k_n * s - k_p*(1-q))

ds/dt = lambda_b * c* P_C *(1-s)- lambda_r *(1-q)*s

dq/dt = (1-q)* k_p * c *(P_C / P_Q)- gamma * q

I have used the ode solver and created two M-files ode.m and lorenz.m

=> Here are my two Matlab M-files. This is my 1st M-file : ode.m which i ran to plot the graph.

  format bank
  close all; 
  clear all; 
  clc; 

  %time interval
  ti=0; 
  tf=140; 
  tspan=[ti tf]; 

  x0=[0.25 0.02 0.98]; %initial vectors

  %time interval of [0 2] with initial condition vector [0.25 0.02 0.02] at time 0.
  options= odeset('RelTol',1e-4, 'AbsTol',[1e-4 1e-4 1e-4]);
  [t,x]= ode23s('lorenz',tspan,x0,options); 

  %Plotting the graphs:
  figure 
  subplot(3,1,1), plot(t,x(:,1),'r'),grid on; 
  title('Lorenz Equations'),ylabel('c'); 

  subplot(3,1,2), plot(t,x(:,2),'b'),grid on; 
  ylabel('s'); 

  subplot(3,1,3), plot(t,x(:,3),'g'),grid on; 
  ylabel('q');xlabel('t')

This is my second M-file which is lorenz.m

  % Creating the MATLAB M-file containing the Lorenz equations.

  function xprime= lorenz(t,x)

   %values of parameters
    I=1200;
    k_f= 6.7*10.^7;
    k_d= 6.03*10.^8; 
    k_n=2.92*10.^9; 
    k_p=4.94*10.^9;
    lambda_b= 0.0087;
    lambda_r =835; 
    gamma =2.74; 
    alpha =1.14437*10.^-3;
    P_C= 3 * 10.^(11);
    P_Q= 2.87 * 10.^(10);    

 % initial conditions
  c=x(1);
  s=x(2);
  q=x(3);

  %Non-linear differential equations.
  % dc/dt= alpha*I*(1-c) + c*(- k_f - k_d - k_n * s - k_p*(1-q))
  % ds/dt = lambda_b * c* P_C *(1-s)- lambda_r *(1-q)*s
  % dq/dt = (1-q)* k_p * c *(P_C / P_Q)- gamma * q

  xprime=[ alpha*I*(1-c) + c*(- k_f - k_d - k_n * s - k_p*(1-q)); lambda_b *(1-s)* c* P_C  - lambda_r *(1-q)*s; (1-q)*k_p * c *(P_C / P_Q)- gamma * q];

Please help me, both M-files codes are working but i want to use function handle (@lorenz) in lorenz.m file because Lorenz isn’t very descriptive of this problem. And also, when i run ode.m file , the values of plot are really small but when i run the lorenz.m file , the values of c,s,q are really big.I want to get values of s and q somewhere between 0 to 1. And value of c should be really big number something 3.5 X10^11. I don't know what is going on?

Manjushree
  • 25
  • 1
  • 6

1 Answers1

1

Your function is incorrect as far as I can see. This line:

xprime=[ alpha*I*(1-c) + c*(- k_f - k_d - k_n * s - k_p*(1-q)); lambda_b * c* P_C - lambda_r *(1-q)*s;   k_p * c *(P_C / P_Q)- gamma * q];

should be:

xprime=[ alpha*I*(1-x(1)) + x(1)*(- k_f - k_d - k_n * x(2) - k_p*(1-x(3))); lambda_b * x(1)* P_C - lambda_r *(1-x(3))*x(2);   k_p * x(1) *(P_C / P_Q)- gamma * x(3)];

You can then get rid of these lines in the function (initial conditions are passed via the call to ode15s):

%initial vectors
  c=0.25; 
  s=0.02;
  q=0.02;
am304
  • 13,758
  • 2
  • 22
  • 40
  • Thank you for your answer. I have changed it but i am still not getting the satisfactory answer. I want to get values of s and q somewhere between 0 to 1. And value of c should be really big number something 3.5 X10^11. – Manjushree Jun 20 '14 at 14:09
  • Check the equations for ds/dt and dq/dt, I think they're wrong (in that the code doesn't correspond to what you have written in your question). You need to start doing your own debugging. – am304 Jun 20 '14 at 14:30
  • i think they are alrite. but how to do debugging in matlab? – Manjushree Jun 20 '14 at 14:39
  • Well, `ds/dt = lambda_b * c* P_C *(1-s)- lambda_r *(1-q)*s` is not the same as `lambda_b * x(1)* P_C - lambda_r *(1-x(3))*x(2)`. The `(1-s)` term is missing. Likewise for dq/dt. – am304 Jun 20 '14 at 14:53
  • xprime=[ alpha*I*(1-c) + c*(- k_f - k_d - k_n * s - k_p*(1-q)); lambda_b * (1-s)* c* P_C - lambda_r *(1-q)*s; (1-q)*k_p * c *(P_C / P_Q)- gamma * q]; This is latest one and i don't think i am missing (1-s) for ds/dt. – Manjushree Jun 20 '14 at 15:01
  • You need to change `c` to `x(1)`, `s` to `x(2)` and `q` to `x(3)`. If after this, it still doesn't give you the expected answer and the equations are correct and correctly implemented, then there's something else wrong, maybe with the parameters. – am304 Jun 20 '14 at 15:10
  • I change c to x(1), s to x(2) and q to x(3) but still getting same answer. So i concern about this line: function xprime= lorenz(t,x) because Sometime it says undefined variable of t when i ran the loren.m file first before ode.m file. And also i want to use function handle (@lorenz) in lorenz.m file because Lorenz isn’t very descriptive of this problem. – Manjushree Jun 20 '14 at 23:00