Решение нелинейной системы из трех дифференциальных уравнений (уравнений Лоренца) с использованием решения ode: ode23s в Matlab

Я пытаюсь решить 3 дифференциальных уравнения (уравнения Лоренца), используя решатель ode: ode23s в Matlab. Вот 3 уравнения Лоренца:

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)- гамма * q

Я использовал решатель ode и создал два M-файла ode.m и lorenz.m.

=> Вот мои два M-файла Matlab. Это мой первый М-файл: ode.m, который я запустил для построения графика.

  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')

Это мой второй М-файл, который называется 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];

Пожалуйста, помогите мне, оба кода М-файлов работают, но я хочу использовать дескриптор функции (@lorenz) в файле lorenz.m, потому что Лоренц не очень описывает эту проблему. А также, когда я запускаю файл ode.m, значения plot действительно малы, но когда я запускаю файл lorenz.m, значения c, s, q действительно велики. Я хочу где-то получить значения s и q от 0 до 1. И значение c должно быть действительно большим числом, что-то 3,5 X10 ^ 11. Я не знаю, что происходит?


matlab nonlinear-optimization matlab-compiler matlab-guide
person Manjushree    schedule 19.06.2014    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Ваша функция неверна, насколько я вижу. Эта строка:

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

должно быть:

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

Затем вы можете избавиться от этих строк в функции (начальные условия передаются через вызов ode15s):

%initial vectors
  c=0.25; 
  s=0.02;
  q=0.02;
person am304    schedule 19.06.2014
comment
Спасибо за ваш ответ. Я изменил его, но я все еще не получаю удовлетворительного ответа. Я хочу получить значения s и q где-то между 0 и 1. И значение c должно быть очень большим числом, что-то 3,5 X10 ^ 11. - person Manjushree; 20.06.2014
comment
Проверьте уравнения для ds/dt и dq/dt, я думаю, что они неверны (поскольку код не соответствует тому, что вы написали в своем вопросе). Вам нужно начать делать свою собственную отладку. - person am304; 20.06.2014
comment
я думаю, что они alrite. а как сделать отладку в матлабе? - person Manjushree; 20.06.2014
comment
Ну, ds/dt = lambda_b * c* P_C *(1-s)- lambda_r *(1-q)*s не то же самое, что lambda_b * x(1)* P_C - lambda_r *(1-x(3))*x(2). Термин (1-s) отсутствует. Аналогично для dq/dt. - person am304; 20.06.2014
comment
xprime=[ alphaI*(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)- гамма * q]; Это последний, и я не думаю, что мне не хватает (1-s) для ds/dt. - person Manjushree; 20.06.2014
comment
Вам нужно изменить c на x(1), s на x(2) и q на x(3). Если после этого он по-прежнему не дает вам ожидаемого ответа, а уравнения верны и правильно реализованы, значит, что-то еще не так, возможно, с параметрами. - person am304; 20.06.2014
comment
Я меняю c на x(1), s на x(2) и q на x(3), но все равно получаю тот же ответ. Поэтому я беспокоюсь об этой строке: function xprime = lorenz (t, x), потому что иногда она говорит о неопределенной переменной t, когда я сначала запускал файл loren.m перед файлом ode.m. А также я хочу использовать дескриптор функции (@lorenz) в файле lorenz.m, потому что Лоренц не очень описывает эту проблему. - person Manjushree; 21.06.2014