Question

Set up and solve a boundary value problem using the shooting method using Matlab

A heated rod with a uniform heat source may be modeled with Poisson equation. The boundary conditions are T(x = 0) = 40 and T(x = 10) = 200 dTf(x) Use the guess values shown below. zg linspace (-200,100,1000); xin-0:0.01:10 a) Solve using the shooting method with f(x) = 25 . Name your final solution "TA" b) Solve using the shooting method with f(x)-0.12x3-2.4x2 + 12x. Name your final solution "TB". Plot both solutions to Poisson's equation so you can compare them and check your boundary conditions.

Answer 1

๖.co) = 40

Matlab code using RK4 for 2nd order differential equation clear all close all %Program for RK4 for question a.. f-e(x,yl,y2) y2; g-e (x,yl,y2) 25; tall guesses for p using shooting method zg=1 in space (-200, 100, 1000) ; for ii-l:length(zg) %function (1) function (ii) x(1,-0 ; y 1 ( 1 )-40;y2(1)-zg ( ii); %initial 加0.01; x_in-x (1); x max-10; n«(x, max-xin) /h ; %number of steps conditions step length %initial x %Final X Runge Kutta 4 iterations for i= 1:n 10-hrg(x(1),y1 ( í ) , y2 (i)); k2-h * f (X ( i) + ( 1 /2 ) *ћ, yl (i ) + ( 1 /2 ) *kly2(i)+(1/2)+11); yl (i+1)-double (yl(i)+(1/6) *(k0+2 k1+2*k2+k3)); y2(ǐ+1)-double (Y2(i) + (1/6) *(10+2+11+2+12+13) ); end yyl(ii)-yl (end)i pinterpl (yyl, zg, 200) clear x; clear yl; clear y2 %solution using value of p x(1,-0 ; y 1 (1)-40;y2 (1)-p; %initial conditions h-0.01; x in-x (1); x max-10; n«(x-max-x-in) /h ; %number of steps step length %initial x %Final x Runge Kutta 4 iterations for i-l:n 10-hrg(x(1),y1 (i),y2(1));

h-0.01 xin-x ( 1 ) ; x, max= 10; n= (x max-x-in) /h; %step length %initial x - %Final x %number of steps iterations %Runge Kutta for i= 1 : n 4 10-htq(x(i),y1 (i),y2 ( i) ) ; k1=h * f ( x ( ī ) + ( 1 /2 ) *h,y1(1)+(1/2)"ko,y2(i)t(1/2)*10); yl (i+1)-double (yl(i)+(1/6) (k0+2*k1+2*k2+k3) y2 (i+1)-double (y2(i)+(1/6) (10+2*11+2 12+13)) end figure (2) %plotting the solution plot (x,yl) title( Plotting of TB(x) using shooting method) xlabel ( ·x' ) ylabel( TB(x) 응욤욤응응욤욤各各욤욤各各욤욤응 End of Code 욤욤各各음음各各음음各各음음各各음

Plotting of TA(x) using shooting method 150 100 150 0 2345 678 910 Plotting of TB(x) using shooting method 200 150 100 50 -50 0 2345 678910

%%Matlab code using RK4 for 2nd order differential equation
clear all
close all
%Program for RK4 for question a..

f=@(x,y1,y2) y2;            %function (i)
g=@(x,y1,y2) 25;            %function (ii)
%all guesses for p using shooting method
zg=linspace(-200,100,1000);
for ii=1:length(zg)
    x(1)=0;y1(1)=40;y2(1)=zg(ii); %initial conditions
    h=0.01;                 %step length
    x_in=x(1);               %Initial x
    x_max=10;           %Final x
    n=(x_max-x_in)/h; %number of steps

    %Runge Kutta 4 iterations
    for i=1:n

        k0=h*f(x(i),y1(i),y2(i));
        l0=h*g(x(i),y1(i),y2(i));
        k1=h*f(x(i)+(1/2)*h,y1(i)+(1/2)*k0,y2(i)+(1/2)*l0);
        l1=h*g(x(i)+(1/2)*h,y1(i)+(1/2)*k0,y2(i)+(1/2)*l0);
        k2=h*f(x(i)+(1/2)*h,y1(i)+(1/2)*k1,y2(i)+(1/2)*l1);
        l2=h*g(x(i)+(1/2)*h,y1(i)+(1/2)*k1,y2(i)+(1/2)*l1);
        k3=h*f(x(i)+h,y1(i)+k2,y2(i)+l2);
        l3=h*g(x(i)+h,y1(i)+k2,y2(i)+l2);
        x(i+1)=x_in+i*h;
        y1(i+1)=double(y1(i)+(1/6)*(k0+2*k1+2*k2+k3));
        y2(i+1)=double(y2(i)+(1/6)*(l0+2*l1+2*l2+l3));
    end
    yy1(ii)=y1(end);
  
end
p = interp1(yy1,zg,200);

clear x; clear y1; clear y2
%solution using value of p

    x(1)=0;y1(1)=40;y2(1)=p; %initial conditions
    h=0.01;                 %step length
    x_in=x(1);               %Initial x
    x_max=10;           %Final x
    n=(x_max-x_in)/h; %number of steps

    %Runge Kutta 4 iterations
    for i=1:n

        k0=h*f(x(i),y1(i),y2(i));
        l0=h*g(x(i),y1(i),y2(i));
        k1=h*f(x(i)+(1/2)*h,y1(i)+(1/2)*k0,y2(i)+(1/2)*l0);
        l1=h*g(x(i)+(1/2)*h,y1(i)+(1/2)*k0,y2(i)+(1/2)*l0);
        k2=h*f(x(i)+(1/2)*h,y1(i)+(1/2)*k1,y2(i)+(1/2)*l1);
        l2=h*g(x(i)+(1/2)*h,y1(i)+(1/2)*k1,y2(i)+(1/2)*l1);
        k3=h*f(x(i)+h,y1(i)+k2,y2(i)+l2);
        l3=h*g(x(i)+h,y1(i)+k2,y2(i)+l2);
        x(i+1)=x_in+i*h;
        y1(i+1)=double(y1(i)+(1/6)*(k0+2*k1+2*k2+k3));
        y2(i+1)=double(y2(i)+(1/6)*(l0+2*l1+2*l2+l3));
    end
    figure(1)
    %plotting the solution
    plot(x,y1)
    title('Plotting of TA(x) using shooting method')
    xlabel('x')
    ylabel('TA(x)')
  
clear x; clear y1; clear y2;

%Program for RK4 for question b..

f=@(x,y1,y2) y2;            %function (i)
g=@(x,y1,y2) 0.12.*x.^3-2.4.*x.^2+12.*x;       
%all guesses for p using shooting method
zg=linspace(-200,100,1000);
for ii=1:length(zg)
    x(1)=0;y1(1)=40;y2(1)=zg(ii); %initial conditions
    h=0.01;                 %step length
    x_in=x(1);               %Initial x
    x_max=10;           %Final x
    n=(x_max-x_in)/h; %number of steps

    %Runge Kutta 4 iterations
    for i=1:n

        k0=h*f(x(i),y1(i),y2(i));
        l0=h*g(x(i),y1(i),y2(i));
        k1=h*f(x(i)+(1/2)*h,y1(i)+(1/2)*k0,y2(i)+(1/2)*l0);
        l1=h*g(x(i)+(1/2)*h,y1(i)+(1/2)*k0,y2(i)+(1/2)*l0);
        k2=h*f(x(i)+(1/2)*h,y1(i)+(1/2)*k1,y2(i)+(1/2)*l1);
        l2=h*g(x(i)+(1/2)*h,y1(i)+(1/2)*k1,y2(i)+(1/2)*l1);
        k3=h*f(x(i)+h,y1(i)+k2,y2(i)+l2);
        l3=h*g(x(i)+h,y1(i)+k2,y2(i)+l2);
        x(i+1)=x_in+i*h;
        y1(i+1)=double(y1(i)+(1/6)*(k0+2*k1+2*k2+k3));
        y2(i+1)=double(y2(i)+(1/6)*(l0+2*l1+2*l2+l3));
    end
    yy1(ii)=y1(end);
  
end
p = interp1(yy1,zg,200);

clear x; clear y1; clear y2
%solution using value of p

    x(1)=0;y1(1)=40;y2(1)=p; %initial conditions
    h=0.01;                 %step length
    x_in=x(1);               %Initial x
    x_max=10;           %Final x
    n=(x_max-x_in)/h; %number of steps

    %Runge Kutta 4 iterations
    for i=1:n

        k0=h*f(x(i),y1(i),y2(i));
        l0=h*g(x(i),y1(i),y2(i));
        k1=h*f(x(i)+(1/2)*h,y1(i)+(1/2)*k0,y2(i)+(1/2)*l0);
        l1=h*g(x(i)+(1/2)*h,y1(i)+(1/2)*k0,y2(i)+(1/2)*l0);
        k2=h*f(x(i)+(1/2)*h,y1(i)+(1/2)*k1,y2(i)+(1/2)*l1);
        l2=h*g(x(i)+(1/2)*h,y1(i)+(1/2)*k1,y2(i)+(1/2)*l1);
        k3=h*f(x(i)+h,y1(i)+k2,y2(i)+l2);
        l3=h*g(x(i)+h,y1(i)+k2,y2(i)+l2);
        x(i+1)=x_in+i*h;
        y1(i+1)=double(y1(i)+(1/6)*(k0+2*k1+2*k2+k3));
        y2(i+1)=double(y2(i)+(1/6)*(l0+2*l1+2*l2+l3));
    end
    figure(2)
    %plotting the solution
    plot(x,y1)
    title('Plotting of TB(x) using shooting method')
    xlabel('x')
    ylabel('TB(x)')

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