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
I'm sorry you are not generating any responses at the moment. Is there any additional information you can share with us? Any new findings?