Location: 12 L Platform 1 model codes @ 4c1ab73f48d7 / USMC / DM+Mijaailovich_huxley_copy / mis.m

Author:
aram148 <a.rampadarath@auckland.ac.nz>
Date:
2021-11-04 16:02:42+13:00
Desc:
Updated USMC-Bursztyn model
Permanent Source URI:
https://staging.physiomeproject.org/workspace/6b0/rawfile/4c1ab73f48d7150c2094cf8017425482a1594ccf/USMC/DM+Mijaailovich_huxley_copy/mis.m
Compare... Download Source 

function [F]=mis(r,N,M,f)
%N number of time points
%trange specifies the interval of time steps
%T time grid
%X0 space grid
T=linspace(0,1000,N);
%X=linspace(xrange(2),xrange(1),M);
X0=linspace(0,1000,M);
%dt=timestep


%Initial conditions
%m10=1;
%m11=0.5;
%m12=2;
%m20=ones(M,1);
%m21=zeros(M,1);
%m22=ones(M,1);
%C=ones(M,1);

R=[r(1);r(2);r(3)];

        
        
        
        
      %  n1=(M10./(q1.*(2*pi).^(0.5))).* exp((-(X-p1).^(2))./(2*(q1.^(2))));
      %loop over number of timesteps to update then run euler
       
        
        %L=lengthvar(t,0.33);
 
    fp1=0.06;
    g1=0.10;
    g2=0.0040;
    g3=0.01;

%gp2=4.*(fp1+gp1);
%gp3=3.*gp1;

        
       r(1)=R(1);
        r(2)=R(2);
        r(3)=R(3);
        %M20=R( (3*M+1):(4*M));
        %M21=R( (4*M+1): (5*M));
        %M22=R( (5*M+1): (6*M));
        %C=R( (6*M+1):  (7*M));
        
      %Mean for cdf
        p1=r(2)/r(1);
      
        %p2=M21(i)/M20(i);
      %Standard deviation for cdf  
        q1=sqrt((r(3)/r(1))-((r(2)/r(1)).^(2)));
       % q2=sqrt((M22(i)/M20(i))-(M21(i)/M20(i)).^(2));
      
    %r,phi and I values for 1st PDE M1_lambda
       r0=-p1./q1;
      % r0=0;
      r1=(1-p1)./q1;
     %r1=1;
      rinf=(Inf-p1)./q1;
      %rinf=Inf;
     phi0=cdf('Normal',r0,p1,q1);
     %phi0=erf(-p1./q1); 
     phi1=cdf('Normal',r1,p1,q1);
     
     %phi1=erf((1-p1)./q1);
     
     phinf=cdf('Normal',rinf,p1,q1);
     
     %phinf=erf((Inf-p1)./q1);
      
      I0=-(exp(-((-p1./q1).^(2))/2))/(sqrt(2*pi));
      I1=-(exp(-((((1-p1)./q1)).^(2))/2))/(sqrt(2*pi));
      I2=-(exp(-(((Inf-p1)./q1)).^(2))/2)/(sqrt(2*pi));
    
      
      
    %Functions for the rhs of the first PDE M1_lambda  
     J0=phi0;
      J01=phi1;
      J0inf=phinf;
      
      J10=((p1.*phi0)+(q1.*I0));
      J11=(p1.*phi1)+(q1.*I1);
      J12=(p1.*phinf)+(q1.*I2);
      
      J20=((p1.^(2)).*phi0)+((2*p1.*q1).*I0)+((q1.^(2)).*(phi0+(r0*I0)));
      J21=((p1.^(2)).*phi1)+((2*p1.*q1).*I1)+((q1.^(2)).*(phi1+(r1*I1)));
      J22=((p1.^(2)))+((2*p1.*q1).*I2)+((q1.^(2)));
     
      J30=(p1.^(3).*phi0)+((3.*p1.^(2).*q1).*I0)+((3.*p1.*q1.^(2)).*((phi0)+(r0*I0)))+((q1.^(3).*(2+r0.^(2)).*I0));
      J31=(p1.^(3).*phi1)+((3.*p1.^(2).*q1).*I1)+((3.*p1.*q1.^(2)).*(phi1+(r1*I1)))+(q1.^(3).*(2+(r1.^(2))).*I1);
      J32=(p1.^(3))+((3.*p1.^(2).*q1).*I2)+((3.*p1.*q1.^(2)));
    
     %Functions defined for the RHS of the second PDE M2_lambda 
      
      %Components for the matrix F that will represent each moment,
      %M1_lambda and M2_lambda
    
%     
     b0=fp1/2;
     b1=fp1/3;
     b2=fp1/4;
     F0=(g2*J0).*r(1)+((fp1+g1)*(J11-J10).*r(1))+g1*(p1-J11).*r(1)+g3*(p1-J11-1+J01).*r(1);
     F1=(g2*J10).*r(1)+((fp1+g1)*(J21-J20).*r(1))+g1*(p1.^(2)+q1.^(2)-J21).*r(1)+g3*(p1.^(2)+q1.^(2)-J21-p1+J11).*r(1);
     F2=(g2*J20).*r(1)+((fp1+g1)*(J31-J30).*r(1))+g1*(p1.^(3)+3*p1*q1.^(2)-J31).*r(1)+g3*(p1.^(3)+3*p1*q1.^(2)-J31-(p1.^(2)+q1.^(2))+J21).*r(1);
      
      % Putting together the matrix F for the RHS of the system of
      % equations
      
      F=[b0-F0;b1-F1;b2-F2];