Question

Discretization, ODE solving, condition number. Consider the differential equation 5y"(x) - 2y'(x) +10y(x)0 on the interval x E [0,10] with boundary conditions y(0)2 and y (10) 3 we set up a finite difference scheme as follows. Divide [0,10] into N-10 sub-intervals, i.e. {xo, X1, [0,1,. 10. Denote xi Xo + ih (here, h- 1) and yi E y(x). Approximate the derivatives as follows X10- 2h we have the following equations representing the ODE at each point Xi ,i = 1, 9 10 y0 2h we omitted the endpoints Xo, x10 since the equations there would require additional points x-1, x11-The boundary conditions y(0)2,y (10) 3 provide two additional equations: yo--2; y1o 3. Adding these and rearranging the general equation we obtain a linear system of 11 equations for the 11 values yo.y,y1o: o-2 a. write the linear system above as Ay-b where У-b0.yı, y10)". (Suggestion: either substitute h 1 or denote the coefficients in the middle equations by p- b. Compute the rank of A. If the matrix A is nonsingular, compute the condition number C. If the system is consistent, solve for y d. Compare yi y(x) obtained this way to the exact solution given below, by computing the Petc) norm of Y as well as that of the error Ywhere Yyexactx) 7x B sin exa) exAc sin(14)

Answer 1

Matlab code for solving Finite difference method ode clear all close all all initial conditions x_in-0;x end-10; y in--2; y end-3; tall step size step length step size h- (x end-x_in)/ (nn+1); all x values x-x in:h:x end; treating tridiagonal matrix for i= 1:nn Creating A matrix elseif i==nn else A(i,it1)- end treating b matrix it i-1 elseif i--nn else bb (1)-0-y in*a; bb (i)-0-y end*ci end end fprint f ( -\tFor fprintf( 'The A disp(A) h-% f\n',h) matrix can be written A=\n') as fprintf(tThe b vector can be written as b-In) disp (bb' ) fprint f ("Rank of A fprintf( 'condition is %f\n", rank (A)) number of A %f\n',cond(A)) Changing the h value

fprintfnltBy changing the h value to 0.02') clear all %all initial conditions x_in-ox_end-10; Yin-2; y_end-3 all step size %step length nn-99; step size h-(x end-x_in)/(nn+1); all x values x®s_in : h : x-end ; %creating tridiagonal matrix for i-1: nn Creating A matrix elseif i-nn else end creating b matrix if i--1 elseif İ_nn else bb( 1 )-0-y-in*a; end end fprintf (n\tFor h-f\n',h) fprint f ("Rank of is %f\n.rank (A)) A fprint f ( 'condition number All y values numerical solutions of %f\n.cond(A)) A y_num(1)-y_in; for i= 1: length(yy) y num(1, i+1)-yy(i); end y_num(1, nn+2)-y_ end;

%Exact solution A=-2 ; B-2*cot ( 14 )+(3*exp (-2))/sin(14); y_ext-(x) exp(x/5)((Acos(7*x/5))+B*sin(7*x/5) %loop for x values for -1:length(x) yy-ext ( 1,1)-yext ( X(1)); end hold on plot (x,y_num, plot (x,yy ext) xlabel(x) ylabel(y(x) titlex vs. y(x) plot' legendNumerical solution', 'Exact solution) %error in exact and numerical solution for h-0.02 error-norm (y_ num-yy_ext); fprint f ('error in exact and numerical solution for h %f is %f. ' r- n',h,error) 용욤욤응용욤욤응용욤욤응용욤 End of Code 各各음음各욤욤음各各음음各욤음음 For h=1.000000 The A matrix can be written as A The b vector can be written as b- 12 Rank of A is 8.000000 Condition number of A Inf By changing the h value to 0.02 For h-o.100000 Rank of A is 99.000000

of A 2445.317704t in tor h-o. 100000 is error in exact and numerical solution for h=0.100000 is 0.281373. x vs. y(x) plot 10 -Numerical solution Exact solution 10 15 0 1 345678 910 Published with MATLABE R2018a

%%Matlab code for solving Finite difference method ode
clear all
close all

%all initial conditions
x_in=0;   x_end=10;
y_in=-2; y_end=3;
%all step size
    %step length
    nn=9;
    %step size
    h=(x_end-x_in)/(nn+1);
    %all x values
    x=x_in:h:x_end;

    %creating tridiagonal matrix
    for i=1:nn
        a=((5/h^2)+(1/h)); b=((-10/h^2)+10); c=((5/h^2)-(1/h));
        %Creating A matrix
        if i==1
            A(i,i)=b;
            A(i,i+1)=c;
        elseif i==nn
            A(i,i)=b;
            A(i,i-1)=a;
        else
            A(i,i)=b;
            A(i,i+1)=c;
            A(i,i-1)=a;
        end

        %creating b matrix
        if i==1
            bb(1)=0-y_in*a;
        elseif i==nn
            bb(i)=0-y_end*c;
        else
            bb(i)=0;
        end
    end
    fprintf('\tFor h=%f\n',h)
    fprintf('The A matrix can be written as A=\n')
    disp(A)
  
    fprintf('\tThe b vector can be written as b=\n')
    disp(bb')
  
    fprintf('Rank of A is %f\n',rank(A))
    fprintf('Condition number of A %f\n',cond(A))
  
  
%Changing the h value
    fprintf('\n\tBy changing the h value to 0.02')
    clear all
    %all initial conditions
    x_in=0;   x_end=10;
    y_in=-2; y_end=3;
    %all step size
    %step length
    nn=99;
    %step size
    h=(x_end-x_in)/(nn+1);
    %all x values
    x=x_in:h:x_end;

    %creating tridiagonal matrix
    for i=1:nn
        a=((5/h^2)+(1/h)); b=((-10/h^2)+10); c=((5/h^2)-(1/h));
        %Creating A matrix
        if i==1
            A(i,i)=b;
            A(i,i+1)=c;
        elseif i==nn
            A(i,i)=b;
            A(i,i-1)=a;
        else
            A(i,i)=b;
            A(i,i+1)=c;
            A(i,i-1)=a;
        end

        %creating b matrix
        if i==1
            bb(1)=0-y_in*a;
        elseif i==nn
            bb(i)=0-y_end*c;
        else
            bb(i)=0;
        end
    end
    fprintf('\n\tFor h=%f\n',h)
  
    fprintf('Rank of A is %f\n',rank(A))
    fprintf('Condition number of A %f\n',cond(A))
  
    %All y values
    yy=A\bb';
    %numerical solutions
    y_num(1)=y_in;
  
    for i=1:length(yy)
        y_num(1,i+1)=yy(i);
    end
  
    y_num(1,nn+2)=y_end;
  
  
    %Exact solution
    A=-2; B=2*cot(14)+(3*exp(-2))/sin(14);
    y_ext=@(x) exp(x/5)*((A*cos(7*x/5))+B*sin(7*x/5));
  
    %loop for x values
    for i=1:length(x)
        yy_ext(1,i)=y_ext(x(i));
    end
  
    hold on
    plot(x,y_num,'r--')
    plot(x,yy_ext)
    xlabel('x')
    ylabel('y(x)')
    title('x vs. y(x) plot')
    legend('Numerical solution','Exact solution')
    %error in exact and numerical solution for h=0.02
    error=norm(y_num-yy_ext);
    fprintf('error in exact and numerical solution for h=%f is %f.\n',h,error)
  
  
    %%%%%%%%%%%%%% End of Code %%%%%%%%%%%%%%%%