非线性优化2 多项式插值，Adam-Bashworth and Euler method

Introduction:

In this lab, I had done some practices about basic structures and operations of MATLAB. Then I learnt some Notation and Repetition of Necessary Mathematical Concepts Notation used throughout the course like:

Vectors and matrices

Linear systems

Derivatives; gradient; Hessian; Jacobian

Taylor expansion

Approximation of derivatives using finite differences.

By using these concepts, I solve some problems in this lab.

Solutions and Analysis:

Lab2.1:

In this problem, we can solve it by using some basic operations

function OrthonormalVectors = GramSchmidt(varargin)     % Every column is an input vector     InputVectors = zeros(2,nargin);     for i=1:nargin        InputVectors(1,i)=varargin{i}(1);     InputVectors(2,i)=varargin{i}(2);     end     [row,col]= size(InputVectors);     OrthonormalVectors = ones(row,col);     OrthonormalVectors(:,1)=InputVectors(:,1);     for i = 2 : col         for j = 1: i-1             InputVectors(:,i)= InputVectors(:,i) - ((OrthonormalVectors(:,j)' * InputVectors(:,i))/(OrthonormalVectors(:,j)' * OrthonormalVectors(:,j))) * OrthonormalVectors(:,j);         end         OrthonormalVectors(:,i)=InputVectors(:,i);     end     % Normalizetion     for i = 1: col         length=norm(OrthonormalVectors(:,i));         for j = 1: row             OrthonormalVectors(j,i)= OrthonormalVectors(j,i)/ length;         end     end   end

 

Then we can use an example to verify this function. By obverting the result of the function. We can find that the function is right.

Lab2.2:

The problem is also very simple; we can complete the answer.

clear all f = @(x,y) x.^2 + y.^2; x0 = 0.5; y0 = -0.5; [xx,yy] = meshgrid(-2:0.1:2); [fx,fy] = gradient(f(xx,yy),0.1); t = (xx == x0) & (yy == y0); at = find(t); fx0 = fx(at); fy0 = fy(at); z = @(x,y) f(x0,y0) + fx0*(x-x0) + fy0*(y-y0); surf(yy,xx,f(xx,yy)); hold on surf(yy,xx,z(xx,yy)); hold on plot3(y0,x0,f(x0,y0)); title('Tangent plane'); xlabel('Y axis'); ylabel('X axis'); zlabel('Function value');

 

 

Lab2.3:

This problem, we need to do some mathematical operations to get some necessary information.

function Polynomial=FitPolynomial(varargin)     InputSetsX = zeros(1,nargin);     InputSetsY = zeros(1,nargin);     for i=1:nargin        InputSetsX(1,i)=varargin{i}(1);     InputSetsY(1,i)=varargin{i}(2);     end     Polynomial=polyfit(InputSetsX,InputSetsY,7); end clear all x1=[0,2,2.8,4,5,6,7]; y1=[0,-1,5,2,-1,5,8]; constant=(FitPolynomial([0,0],[2,-1],[2.8,5],[4,2],[5,-1],[6,5],[7,8])); [~,col]=size(constant); syms x ; Polynomial= (constant(1)*x.^7.+constant(2)*x.^6.+constant(3)*x.^5.+constant(4)*x.^4.+constant(5)*x.^3.+constant(6)*x.^2.+constant(7)*x.^1.+constant(8)); x=0:0.001:7; dPoly=diff(Polynomial,1); fun=@(x)x.*4.284434610872131e+1-x.^3.*1.172790362641904e+1+x.^4.*3.976413349884542-x.^5.*4.829189890543072e-1+x.^6.*1.965463551562472e-2-3.0217193158606e+1; [points(1),values(1)]= fzero(fun,[0.2,2]); [points(2),values(2)]= fzero(fun,[2,4]); [points(3),values(3)]= fzero(fun,[4,5.3]); [points(4),values(4)]= fzero(fun,[5.3,6.8]); dPoly=eval(dPoly); Polynomial = eval(Polynomial); plot(x,dPoly,'b','LineWidth',2); title('Polynomial Interpolation'); hold on plot(x1,y1,'bo','LineWidth',2); hold on plot(x,Polynomial,'k','LineWidth',2); hold on plot([0,7],[0,0],'r','LineWidth',2); hold on for i=1:4     values(i)=polyval(constant,points(i)); end plot(points,values,'rx','LineWidth',2); hold on for i=1:2     plot([points(i*2-1),points(i*2-1)],[values(i*2-1),0],'g--','LineWidth',2);     hold on     plot([points(i*2),points(i*2)],[0,values(2*i)],'g--','LineWidth',2);     hold on end legend('levelFirst-order derivative of p','Interpolation polynomial p','Data points','Zero Level','Locations of local maxima/minima');

Lab2.4:

In this problem, we need to understand the principle of two kinds of methods to solve ODE. The first one is the Euler’s method which is the easiest one. The other one is Adam-Bashworth which is a little bit difficult but has better performance on accuracy.

clear all syms x y % Range of x xspan=[0,5]; % Step size h=0.1; % Initial values y0=1, x0=0; y(1)=y0;x(1)=x0; % No. of steps steps=(xspan(2)-xspan(1))/h; %% Define f(x,y) %% f= @(x,y)y; %% Euler's Method   for j=2:steps+1         x(j,1)=x(j-1)+h;         y(j,1)=y(j-1,1)+h*f(x(j-1,1),y(j-1,1));     end euler=vpa([x y],5); %% Adam-Bashworth for k=5:steps+1 x(k,1)=x(k-1)+h;    y(k,1)=y(k-1) +(h/24)*( -9*f(x(k-4),y(k-4)) +37*f(x(k-3),y(k-3))...                         -59*f(x(k-2),y(k-2)) +55*f(x(k-1),y(k-1))); end p=y; Adam_Bashworth=vpa([x p],5);   syms s t s=exp(t); t=linspace(0,5,1000); s=eval(s); figure plot(t,s,'-b','LineWidth',1)   %plot(Exact_Sol(:,1),Exact_Sol(:,2),'b'); hold on plot(euler(:,1),euler(:,2),'k--'); hold on plot(Adam_Bashworth(:,1),Adam_Bashworth(:,2),'r'); legend('Exact solution','Numerical solution (Euler)','Numerical solution (Adams-Bashforth)');

The result of this function is shown like following:

Figure 1 Case1 with h=0.2

Figure 2 Case1 with h=0.1

Figure 3 Case1 with h=0.05

Figure 4 Case1 with h=0.01

 

Figure 5 Case2 with h=0.2

Figure 6  Case2 with h=0.1

Figure 7 Case2 with h=0.05

Figure 8  Case2 with h=0.01

By obverting the result of two kinds of methods.

We can find that for the performance:

The Adam_Bashworth had better effect.

And when we use the same method, the performance is better with the interval get smaller.