Question

Problem 2: consider the non-dimensionalized model for dynamic protein concentration: 2.g(0)-c A) Use qualitative analysis (plot the r.h.s. numerically) and find the steady states for this equation for a 0.2, b 0.01. Identify which steady states are stable and which are not. B) Solve the equation for these parameters numerically using the initial condition c 10. Solve it one more time for c 2. (Choose the interval on which you solve the ODE to be able to explain the result.) Plot the results, and explain them in 2-3 sentences.

Answer 1

dg/dt = 1 - g + ag^2/1 + bg^2, g(0) = c (A) here a = 0.2, b = 0.01 therefore dg/dt = 1 - g + 0.2 times g^2/1 + 0.01g^2 plotting of the functions f(g) = 1 - g + 0.2 times g^2/1 + 0.01g^2 here the roots are A = 1.36675 B = 4.99 C = 14.63 here g = 1.36675 for point A is stable as f(A - A) > 0 f(A + A) < 0 for g = 4.99 for point B is unstable as f(B - A) < 0 f(B + A) > 0 and g = 14.63 for point c is stable as f(C - A) > 0 f(C + A) < 0 here C = 10 it will go to stable point 14.63 and converges for C = 2 it will go to stable point 1.3667 and converges \r\n %%Matlab code for Direction field and solution of ODE clear all close all %function of 1st order differential eqn

%%Matlab code for Direction field and solution of ODE
clear all
close all
%function of 1st order diffrential eqn
f = @(g) 1-g+((0.2*g.^2)./(1+0.01.*g.^2));
%plotting of R.H.S of the function
gg=-5:0.1:20;
f_gg=f(gg);
plot(gg,f_gg)
grid on
title('plotting of R.H.S of the function')
xlabel('g')
ylabel('f(g)')

%Finding the location for which function f crosses zero
%Here we will going to use Bisection method to find the root
%all initial guess
a=[0 7 7];
b=[3 3 20];
fprintf('For the function \n')
disp(f)
fprintf('The roots are \n')
fprintf('\ta\tb\troot\n')
for i=1:length(a)
    [root1,count1]=bisection_method(f,a(i),b(i));
    fprintf('\t%0.2f\t%0.2f\t%f\t\n',a(i),b(i),root1)
end
fprintf('Here a and b are intervals\n')

%function of 1st order diffrential eqn
fun = @(t,g) 1-g+((0.2*g.^2)./(1+0.01.*g.^2));
%solution of differential equation using Runge Kutta 4
%step size
    h=0.01;
%all final t steps
    t1=0;tn=20;
    t_in=t1;     %Initial t
    t_max=tn;    %Final t
    %Runge Kutta 4 iterations
    n=(t_max-t_in)/h;
    g_rk(1)=10; %initial g
    t_rk(1)=t1;
    %RK4 iterations
    for i=1:n
     
        k0=h*fun(t_rk(i),g_rk(i));
        k1=h*fun(t_rk(i)+(1/2)*h,g_rk(i)+(1/2)*k0);
        k2=h*fun(t_rk(i)+(1/2)*h,g_rk(i)+(1/2)*k1);
        k3=h*fun(t_rk(i)+h,g_rk(i)+k2);
        t_rk(i+1)=t_in+i*h;
        g_rk(i+1)=double(g_rk(i)+(1/6)*(k0+2*k1+2*k2+k3));
      
    end
figure(2)
plot(t_rk,g_rk)
ylabel('g(t)')
xlabel('t')
title('solution plot for differential equation c=10')
grid on
clear t_rk; clear g_rk
%solution of differential equation using Runge Kutta 4
%step size
    h=0.01;
%all final t steps
    t1=0;tn=20;
    t_in=t1;     %Initial t
    t_max=tn;    %Final t
    %Runge Kutta 4 iterations
    n=(t_max-t_in)/h;
    g_rk(1)=2; %initial g
    t_rk(1)=t1;
    %RK4 iterations
    for i=1:n
     
        k0=h*fun(t_rk(i),g_rk(i));
        k1=h*fun(t_rk(i)+(1/2)*h,g_rk(i)+(1/2)*k0);
        k2=h*fun(t_rk(i)+(1/2)*h,g_rk(i)+(1/2)*k1);
        k3=h*fun(t_rk(i)+h,g_rk(i)+k2);
        t_rk(i+1)=t_in+i*h;
        g_rk(i+1)=double(g_rk(i)+(1/6)*(k0+2*k1+2*k2+k3));
      
    end
figure(3)
plot(t_rk,g_rk)
ylabel('g(t)')
xlabel('t')
title('solution plot for differential equation c=2')
grid on
%code for direction field
figure(4)
hold on
%plotting dir field for differential equation
dirfield(fun,0:1:20,-2:1:20)

%Findind solution for all initial condition -2:2
for y0=-2:1:20
%solution using matlab inbuilt ode45 function
[xs,ys] = ode45(fun,[0,20],y0);
%plotting of solution
   plot(xs,ys,'b')
end
hold off
xlabel('t')
ylabel('g(t)')
title('Direction field and g(t) vs t plot')

%%Matlab function for bisection method
function [root,count]=bisection_method(fun,a,b)
%f(a1) should be positive
%f(b1) should be negative
count=0;
k=10;
%loop for bisection method
while k>10^-6
    count=count+1;
    xx(count)=(a+b)/2;
    mm=double(fun(xx(count)));
    if mm>=0
        a=xx(count);
    else
        b=xx(count);
    end
    k=abs(a-b);
end
root=xx(end);
end

%function for plotting direction field
function dirfield(f,tval,yval)
% dirfield(f, t1:dt:t2, y1:dy:y2)
%
%   plot direction field for first order ODE y' = f(t,y)
%   using t-values from t1 to t2 with spacing of dt
%   using y-values from y1 to t2 with spacing of dy
%
%   f is an @ function, or an inline function,
%     or the name of an m-file with quotes.
%
% Example: y' = -y^2 + t
%   Show direction field for t in [-1,3], y in [-2,2], use
%   spacing of .2 for both t and y:
%
%   f = @(t,y) -y^2+t
%   dirfield(f, -1:.2:3, -2:.2:2)

[tm,ym]=meshgrid(tval,yval);
dt = tval(2) - tval(1);
dy = yval(2) - yval(1);
fv = vectorize(f);
if isa(f,'function_handle')
fv = eval(fv);
end
yp=feval(fv,tm,ym);
s = 1./max(1/dt,abs(yp)./dy)*0.35;
h = ishold;
quiver(tval,yval,s,s.*yp,0,'.r'); hold on;
quiver(tval,yval,-s,-s.*yp,0,'.r');
if h
hold on
else
hold off
end
axis([tval(1)-dt/2,tval(end)+dt/2,yval(1)-dy/2,yval(end)+dy/2])
end

%%%%%%%%%%%%%%%%%%% End of Code %%%%%%%%%%%%%%%%