"find $p(1)$, $p(2)$, $p(3)$ and $p(4)$ which fit\r\n$N = [9.045015676e6 1.899453292e7 1.989903449e7 5.427009406e7 ...\r\n 9.045015676e7 2.532604389e8 3.979806897e8 4.341607524e8 ...\r\n4.432057681e8 5.246109092e8 6.150610660e8];$\r\nby use fminsearch code in Matlap.\r\n\\begin{equations}\r\ndydt(1) = 0.515*y(1)*(1-y(1)*1e-9)/(1 + 8.164*y(1))- p(1)*y(3)*y(1)- p(2)*y(4)*y(1)- p(3)*y(5)*y(1)+ p(4)*y(6)*y(1) - 0.01525*y(2); \r\n\r\ndydt(2) = -0.69*y(2) + 10*(heaviside(t-10) - heaviside(t-10.2)); \r\ndydt(3) =0.09*y(3)*y(5)*(1-(1/100000000)*y(3)-(1/100000000)*y(4))-0.03*y(3)+0.008*y(5)-(1/10000000)*y(3)*y(1)+0.00594*y(3)*y(2);\r\ndydt(4) =0.09*y(4)*y(5)*(1-(1/100000000)*y(3)-(1/100000000)*y(4))-0.03*y(4)+0.008*y(6)-(1/10000000)*y(4)*y(1)+0.00594*y(2)*y(4);\r\ndydt(5) = 0.02*y(5)*(1-(1/1000000000)*y(5)-(1/1000000000)*y(6))-0.07*y(5)+0.001*y(3)+0.09*y(6)-(1/10000000)*y(5)*y(1)+0.01*y(5)*y(2);\r\ndydt(6) = 0.02*y(6)*y(4)*(1-(1/1000000000)*y(5)-(1/1000000000)*y(6))-0.11*y(6)+0.001*y(4)+0.05*y(5)+(1/10000000)*y(6)*y(1)-0.01*y(6)*y(2); \r\n\\end{equations}\r\n\r\nSolution to a System of ODEs May 13, 2017 The solution to the given system of ordinary di?erential equations depend on four parameters, that is given by the cubic approximationto the data vector N. The coe?cients of the third degree polynomial can be found by the Matlab function fminsearch. This routine ?nds the solution to the unconstrained opti- mization problem of minimizing the squared errorbetween the function and the data. The result is shown in Figure 1. Figure 1: Least squares ?t of polynomial of degree three. By requiring that the all p = 0, a non-linear optimization routine for con- i strained problems (fmincon) is invoked. Again the error of the ?tted curve and the data is minimized, now with the constraint that p = 0. The modi?ed i polynomial is shown in Figure 2. We are given a system of ODEs (in Matlab notation), dydt(1) = 0.515*y(1)*(1-y(1)*1e-9)/(1 + 8.164*y(1)) - p(1)*y(3)*y(1) - p(2)*y(4)*y(1) - p(3)*y(5)*y(1) + p(4)*y(6)*y(1) 1Figure 2: Optimization of parameters, forcing pi= 0. - 0.01525*y(2); dydt(2) = -0.69*y(2) + 10*(heaviside(t-10) - heaviside(t-10.2)); dydt(3) = 0.09*y(3)*y(5)*(1-(1/100000000)*y(3) - (1/100000000)*y(4)) - 0.03*y(3)+0.008*y(5) - (1/10000000)*y(3)*y(1) + 0.00594*y(3)*y(2); dydt(4) = 0.09*y(4)*y(5)*(1-(1/100000000)*y(3) - (1/100000000)*y(4)) - 0.03*y(4)+0.008*y(6) - (1/10000000)*y(4)*y(1) + 0.00594*y(2)*y(4); dydt(5) = 0.02*y(5)*(1-(1/1000000000)*y(5) - (1/1000000000)*y(6)) - 0.07*y(5) + 0.001*y(3) + 0.09*y(6) - (1/10000000)*y(5)*y(1) + 0.01*y(5)*y(2); dydt(6) = 0.02*y(6)*y(4)*(1-(1/1000000000)*y(5) - (1/1000000000)*y(6)) - 0.11*y(6)+0.001*y(4)+0.05*y(5) + (1/10000000)*y(6)*y(1)-0.01*y(6)*y(2); The 6-dimensional system of ODEs is very sensitive to the values of p . i Indeed, using the ?rst approximation in the system, one is faced with severe numerical di?culties, where the solution causes numerical over?ow after only a short time interval. The solution of the system using the modi?ed (optimized) values of p is i 2shown in Figure 3. The solution is computed using the Matlab function ode45 with initial value (50000,0,0,0,100,0) and maximum time step dt = 0.1. To avoid the initial large spikes, the solution is plotted for t= 1. Figure 3: Solution to scaled system of ODEs. The e?ect of the step function at t = 10 is clearly visible. The higher dimensions are small compared to y and y , and tend to zero as t?8. 1 2 1 MATLAB Code % Data file for assignment N = [9.045015676e6 1.899453292e7 1.989903449e7 5.427009406e7 ... 9.045015676e7 2.532604389e8 3.979806897e8 4.341607524e8 ... 4.432057681e8 5.246109092e8 6.150610660e8]; x = [0:length(N)-1]; % Curve fit using cubic approximation and % Script for solving a system of ODEs matlab = true; FourData guess = [1e-9 1e-9 1e-9 1e-9]; disp(’First approximation of p:’) p0 = FourApprox(guess) 3ya = p0(1)+p0(2).*x+p0(3).*x.^2+p0(4).*x.^3; figure; plot(x, ya, x, N, ’o’) title(’Cubic Least Squares Approximation’) xlabel(’x’) ylabel(’Value (N)’) disp(’Optimization with p>=0’) if (matlab) [err, p] = fmincon(@poly3err, guess, [], [], [0 0 0 0], []) else [p, err] = sqp(guess, @poly3err, [], [], [0 0 0 0], []) end pval = p(1)+p(2).*x+p(3).*x.^2+p(4).*x.^3; figure; plot(x, pval, x, N, ’o’) title(’Cubic Least Squares Modification’) xlabel(’x’) ylabel(’Value (N)’) y0 = [50000 0 0 0 100 0]’; dt = 0.1; t0 = [0:dt:20]; options = odeset(’InitialStep’,dt,’MaxStep’,dt); if (matlab) [t,y] = ode45(@FourPe, t0, y0, [], p, options); else [t,y] = ode45(@FourPe, t0, y0, p, options); end startidx = 10; tplot = t(startidx:length(t)); yplot = y(startidx:length(y),:); figure; plot(tplot,yplot) title(’Solutions to System of ODEs’) xlabel(’t’) ylabel(’y’) legend(’y(1)’,’y(2)’,’y(3)’,’y(4)’,’y(5)’,’y(6)’) 4% Function file for cubic polynomial fit to data function p = FourApprox(guess) matlab = true; FourData if (matlab) p = fminsearch(@poly3err, guess); else p = polyfit(x, N, 3); p = fliplr(p); end end % Function file for error of cubic polynomial fit function err = poly3err(p) FourData pval = p(1)+p(2).*x+p(3).*x.^2+p(4).*x.^3; err = 0; for i = 1:length(N) err = err + (pval(i) - N(i))^2; end end % Function file system of ODEs function dydt = FourPe(t,y,p) dydt(1) = 0.515*y(1)*(1-y(1)*1e-9)/(1 + 8.164*y(1)) - p(1)*y(3)*y(1) - p(2)*y(4)*y(1) - p(3)*y(5)*y(1) + p(4)*y(6)*y(1) - 0.01525*y(2); dydt(2) = -0.69*y(2) + 10*(heaviside(t-10) - heaviside(t-10.2)); dydt(3) = 0.09*y(3)*y(5)*(1-(1/100000000)*y(3) - (1/100000000)*y(4)) - 0.03*y(3)+0.008*y(5) - (1/10000000)*y(3)*y(1) + 0.00594*y(3)*y(2); dydt(4) = 0.09*y(4)*y(5)*(1-(1/100000000)*y(3) - (1/100000000)*y(4)) - 0.03*y(4)+0.008*y(6) - (1/10000000)*y(4)*y(1) + 0.00594*y(2)*y(4); dydt(5) = 0.02*y(5)*(1-(1/1000000000)*y(5) - (1/1000000000)*y(6)) - 0.07*y(5) + 0.001*y(3) + 0.09*y(6) - (1/10000000)*y(5)*y(1) + 0.01*y(5)*y(2); dydt(6) = 0.02*y(6)*y(4)*(1-(1/1000000000)*y(5) - (1/1000000000)*y(6)) - 0.11*y(6)+0.001*y(4)+0.05*y(5) + (1/10000000)*y(6)*y(1)-0.01*y(6)*y(2); end 5"