Solving Non-linear System of 3 diff eqns using ode23s in matlab

Join the discussion
Registration is free. Ask a follow-up in this thread, or start your own.
1 reply · 4K views
wel
Gold Member
Messages
36
Reaction score
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.
Code:
      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

Code:
      % 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?
 
Physics news on Phys.org