Shooting Method For Boundary Value Problem Solving

Shooting Method for Boundary Value Problem Solving ဒီတခါEတာ့ စာဖတ္သူမ်ားကUI Boundary Value Problem Solving method တစ္ခUၿဖစ္တဲ့ Shooting method ဆUIတာနဲ႔ မIတ္ဆက္Eပးပါမယ္။ Aရင္က Eရးဘူးတဲ့ Continuous Method နဲ႔ Compare လUပ္ၿပီးA သUံးၿပUလUI႔ ရEAာင္လUI႔ပါ။ matlab မွာEတာ႔ bvp4c, bvp5c ဆUIတဲ့ functions EတြရွIပါတယ္။ ဒါEပမယ့္ AသUံး ၿပUရတာ မလြယ္ကူတဲ့Aၿပင္ ရလာမယ့္ AEၿဖကUIလဲ ဘယ္လUIရလာတယ္ ဆUIတာၿပန္ရွင္းၿပဖUI႔ AEတာ္ခက္ပါလIမ့္ မယ္။ Propblem တစ္ခUကUI solve လUပ္မယ္ဆUIရင္ Aဲဒီ problem ရဲ့ AEၿဖကUI ရEAာင္ ရွာႏUIင္ဖUI႔ Aၿပင္ AဲဒီA Eၿဖဘယ္လUIရလာတယ္ဆUIတာ ကUIသIဖUI႔ကလဲ AEရးႀကီးပါတယ္။ ဒီ Posts EတြEရးရတဲ့ ရည္ရြယ္ခ်က္ကလဲ ODE BVP တစ္ခက U UI နားလည္ၿပီးEတာ့ ဘယ္လUIနည္းEတြနဲ႔ Solve လUပ္လUI႔ ရတယ္ဆUIတာ သIEစခ်င္တဲ့ ရည္ ရြယ္ခ်က္ပါ။ ဒါမွလဲ Engineering Problems Eတြမွာ လUIAပ္ရင္ လUIAပ္သလUI AသUံးၿပU ႏIUင္မွာၿဖစ္ပါတယ္။ ကဲ.. စလUIက္ႀကရEAာင္။ Example

to=0, tf=1, y1o=y2=1, y1f=2, y2f=3; Problem က y1,y2,y3,y4 ဆUIတဲ့ Differential equations EလးခUရွIမယ္။ t initial and t final ရွIမယ္။ y1 and y2 ရဲ့ initial value EတြသIမယ္။ လUIခ်င္တဲ့ y1 and y2 ရဲ့ final value EတြကUI ရဖUI႔ Aတြက္ y3 and y4 ရဲ့ Initial value ကUI ရွာဖUI႔ပါဘဲ။ ကဲ.. Shooting method ကUI သUံးၿပီးရွာႀကည့္ႀကရEAာင္။ Shooting method က ကၽြန္Eတာ္တUI႔ လUIခ်င္တဲ့ initial values ကUI guess လUပ္ၿပီး Eပးထားတဲ့ ODE ကUI integrate လUပ္ကာ the system of nonlinear equations Aသြင္Eၿပာင္းၿပီး Solve လUပ္ၿခင္းဘဲ။ Eပးထားတဲ့ ODE ကUI matlab မွာ Eရးႀကည့္ရEAာင္။ function F = dEqs(t,y) % Differential equations. F = zeros(4,1); F(1) = y(3)^2; F(2) = sqrt(y(4))-y(1); F(3) = y(2)+y(1); f(4) = y(1)-y(3);

ဒီ ODE ကUI integrate လUပ္ဖUI႔ရန္Aတြက္ initial value က y1 and y2 ပဲ သIတယ္ဆUIတဲ့ Aတြက္ y3 and y4 တUI႔ရဲ့ Initial value ကUI u(1) and u(2) လUI႔ assign လUပ္ၿပီး Initial Condition functionကUI Eရးပါမယ္။ function y = inCond(u) % Initial conditions; u(1) y = [1 1 u(1) u(2)]; % and u(2) are unknowns.

Initial condition ရၿပီ ဆUIEတာ့ Eပးထားတဲ့ ODE ကUI Integrate လUပ္ၿပီး the system of nonlinear equations Aသြင္ရဖUI႔ function Eရးႀကည့္ရEAာင္။ function r = residual(u) % Bounday residuals. global XSTART XSTOP r = zeros(length(u),1); x = XSTART; [xSol,ySol] = ode45(@dEqs,[x XSTOP],inCond(u)); [m,n] = size(ySol); r(1) = ySol(m,1)-2; r(2) = ySol(m,2)-3;

ဒီ function မွာမသIတာက u(1) and u(2) ၿဖစ္ၿပီး Eၿပလည္ရမယ့္ equations က r(1) and r(2) ၿဖစ္လာပါၿပီ။ A ဲEတာ့ two unknowns and two equations ပUံစံ ရလာၿပီဆUIEတာ့ ဒါကUI ရွင္းဖUI႔ solving the system of nonlinear equations ကUI ႀကည့္ႀကရEA ာင္။ solving system of nonlinear equations ဆUIတာက F(x)=0 ဆUIတာကUI Eၿဖရွင္းဖUI႔ပါဘဲ။ Matlab မွာEတာ့ ဒါကUIEၿဖရွင္းဖUI႔ fsolve ဆUIတဲ့ function ရွIပါတယ္။ fsolve မွာက built-in A Eနနဲ႔ သUံးထားတာက trust region method ပါ။ သူ႔ရဲ့

built-in algorithm နဲ႔မွ A ဆင္မEၿပဘူး

ဆUIရင္ options မွာ Gauss-newton or Levenberg-Marquardt A သUံးၿပUလUI႔ရEA ာင္ ဝင္ၿပီး change ယူလUI႔ရ ပါတယ္။ ဒီEနရာမွာ စာဖတ္သူမ်ားကUI root finding နဲ႔ ပတ္သက္ၿပီးEတာ့ ပUIၿပီး နားလည္သEဘာEပါက္EစဖUI႔ newton-raphson method ကUIA သUံးၿပUၿပီး ရွာႀကည့္ႀကရEA ာင္။ newton-raphson method ကUIEတာ့မရွင္း Eတာ့ပါဘူးEနာ္။ Advanced Engineering mathematics မွာ A ားးလUံးသIကၽြမ္းရင္းႏွီးခဲ့ၿပီးသားလUI႔ ထင္လUI႔ပါ။ Newton-rapson method ကUI matlab မွာEA ာက္ပါA တUIင္းEရးလUI႔ရပါတယ္။ function root = newtonRaphson2(func,x,tol) % Newton-Raphson method of finding a root of simultaneous % equations fi(x1,x2,...,xn) = 0, i = 1,2,...,n. % USAGE: root = newtonRaphson2(func,x,tol) % INPUT: % func = handle of function that returns[f1,f2,...,fn]. % x = starting solution vector [x1,x2,...,xn]. % tol = error tolerance (default is 1.0e4*eps). % OUTPUT: % root = solution vector. if nargin == 2; tol = 1.0e4*eps; end

if size(x,1) == 1; x = x'; end % x must be column vector for i = 1:150 [jac,f0] = jacobian(func,x); if sqrt(dot(f0,f0)/length(x)) < tol root = x; return end dx = jac\(-f0); %dx=-1e-3; x = x + dx; if sqrt(dot(dx,dx)/length(x)) < tol*max(abs(x),1.0) root = x; return end end error('Too many iterations')

function [jac,f0] = jacobian(func,x) % Returns the Jacobian matrix and f(x). h = 1.0e-4; n = length(x); jac = zeros(n); f0 = feval(func,x) for i =1:n temp = x(i); x(i) = temp + h; f1 = feval(func,x); x(i) = temp; jac(:,i) = (f1 - f0)/h; end

ဒီEနရာမွာ newton-raphson မွာ လUIA ပ္တဲ့ jacobian function ကUI finite difference method ကUI သUံးၿပီး numerical Eရးထားပါတယ္။A ဲ.. Newton-raphson function လဲ ရၿပီဆUIEတာ့ A ထက္မွာEရးထားတဲ့ two unknowns two equations function ကUI ရွင္းႀကည့္ႀကရEA ာင္။ function shoot4nl % Shooting method for nonlinear 4th-order boundary % value problem in Example. global XSTART XSTOP % Make these params. global. XSTART = 0; XSTOP = 1; % Range of integration. u = [-1 1]; % Trial values of u(1)

% and u(2). x = XSTART; u = newtonRaphson2(@residual,u)

ဒီ function ကUI run လUIက္တာနဲ႔ ကၽြန္Eတာ္တUI႔လUIခ်င္တဲ့ y3 and y4 ရဲ့ initial value ကUI EA ာက္ပါA တUIင္း ရမွာၿဖစ္ပါတယ္။ u= -1.3777 11.4519 y3 and y4 ရဲ့ initial value EတြရလာၿပီဆUIEတာ့ ရလာတဲ့ တန္ဖUIးEတြဟာ ကၽြန္Eတာ္တUI႔ ရွင္းခ်င္တဲ့ Boundary Value Problem ကUI Eၿပလည္EစသလားဆUIတာ စစ္ႀကည့္ရEA ာင္။ [xSol,ySol] = ode45(@dEqs,[x XSTOP],inCond(u)); y1=ySol(:,1); y2=ySol(:,2); plot(xSol,y1) grid on figure plot(xSol,y2)

Fig(1) change of y1

U U

fig(2).change of y2

fig(1) and fig(2) ကUI ႀကည့္ၿခင္းA ားၿဖင့္ ကၽြန္Eတာ္တUI႔ ရွင္းခ်င္တဲ့ BVP Eၿပလည္EA ာင္ရွင္းနIUင္Eႀကာင္းသI ႏIUင္ပါတယ္။ ဒီEနရာမွာ Newton-raphson method

မသUံးဘဲ matlab က fsolve function ကUIသUံးၿပီး

ရွာႀကည့္ရEA ာင္။ u တန္ဖUIးဟာ EA ာက္ပါA တUIင္း ရရွIမွာကUIEတြ႔ရပါတယ္။ u= -1.2250 11.0197 ဒီ y3 and y4 ရဲ့ initial တန္ဖUIးEတြသUံးၿပီးEတာ့ ODE ကUI integrate လUပ္ႀကည့္မယ္ဆUIရင္ EA ာက္ပါA တUIင္းရ ရွIမွာၿဖစ္ပါတယ္။

Fig(3).change of y1

Fig(4).change of y4 fig(3) and fig(4) ကUI ႀကည့္ၿခင္းA ားၿဖင့္ ကၽြန္Eတာ္တUI႔ ရွင္းခ်င္တဲ့ BVP Eၿပလည္EA ာင္ရွင္းနIUင္Eႀကာင္းသI ႏIUင္ပါတယ္။ ဒီEနရာမွာ ထပ္Eၿပာခ်င္တာက solving system of nonlinear equations A Eႀကာင္းပါ။ system of nonlinear equations ကUI ရွင္းဖUI႔A တြက္ numerial methods Eတြ A မ်ားႀကီးရွIပါတယ္။ A ဲဒါEတြ ကEတာ့

(1) steepest desent method (2) bisection method (3) newton method (4) newton-raphson method (5) gauss-newton method (6) quasi-newton method

(7) Levenberg-Marquardt method (8) dogleg method စသည္ၿဖင့္Eပါ့Eလ။ ဒီEနရာမွာ စာဖတ္သူA Eနနဲ႔ မIမIရဲ့ problem A Eပၚမူတည္ၿပီး ဘယ္ method ကUI A သUံး ၿပUမလဲဆUIတာကUI ခ်င့္ခ်Iန္ရမွာၿဖစ္ပါတယ္။ matlab A သUံးၿပUသူEတြA Eနနဲ႔ကEတာ fsolve ကUI ပUIင္ပUIင္ႏIUင္ႏUIင္သUံး တတ္ရင္A ဆင္Eၿပမွာပါ။ သူသUံးထားတဲ့ background algorithm Eတြက Eတာ္Eတာ္ကUI strong ၿဖစ္ပါတယ္။ Eတာ္Eတာ္မ်ားမ်ား problems Eတြ convergence ၿဖစ္ပါတယ္။ ဒါEပမယ့္ စာEရးသူ Myanmar Astronomy ထဲမွာ တင္Eပးထားတဲ့ presentation for Shanghai ဆUIတဲ့ထဲက BVP ကUIEတာ့ convergence မၿဖစ္ပါဘူး။ differential equations Eတြက A ရမ္းရွUပ္Eထြးတာရယ္၊ boundary conditions Eတြမ်ားတဲ့ A ခါက်Eတာ့ initial တန္ဖUIးက A ရမ္းကUI sensitive ၿဖစ္လြန္းEနလUI႔ပါ။ ဒါEပမယ့္ A ဒီ Problem ကUI စာEရးသူ A ရင္ post က တင္ထားတဲ့ continuous method သUံးၿပီးEၿဖရွင္းတဲ့A ခါမွာEတာ့ convergence ၿဖစ္ပါတယ္။ ဒါEပမယ့္ core2duo prosessor, RAM 4Gb နဲ႔ ၁၀နာရီEလာက္Eတာ့ Run ယူရပါတယ္။ A ဲဒါက လြန္ခဲ့တဲ့ ၆လEလာက္ တUန္းက solve လUပ္လUI႔ရခဲ့တာEပါ့။ A ခUEတာ့ A ဲလUI ၁၀နာရီEလာက္ run စရာမလUIဘဲ ၁၀မIနစ္Eလာက္ run ရUံနဲ႔ A Eၿဖရမယ့္ method တခUEတာ့ ရပါၿပီ။ A ဲဒီ method Eတာ့ powell’s hybrid method လUI႔Eခၚပါတယ္။ A ခ်Iန္ရမယ္ဆUIရင္Eတာ့ next post Eတြမွာ Eရးၿဖစ္မယ္ထင္ပါတယ္။ ကၽြန္Eတာ္ကEတာ့ bluephoenix လUI A လြမ္းကဗ်ာEတြမEရးတတ္Eတာ့ ၈ဏန္းသခၤ်ာEတြပဲEရးၿဖစ္ မယ္ထင္ပါတယ္။ ငယ္ငယ္တUန္းက ဟာသဝUIင္းလUပ္ခဲ့ႀကတဲ့ ဝါက်Eလး တစ္ခUရွIပါတယ္။ ‘ A လြမ္းEတြကUI ပန္းEတြနဲ႔ သက္Eသၿပရရင္ၿဖင့္ မီး တUI႔Eနတဲ့ ကမာၻႀကီးဟာ ပန္းကမာၻႀကီး ၿဖစ္သြားပါလIမ့္မယ္ ကUIကUI’ ဆUIတာပါ။ A ခUEတာ့

ကၽြန္Eတာ့္A Eနနဲ႔

A ဲ ဒီ

ဝါက်ကUI

ဒီလUIEၿပာင္းခ်င္တယ္။

‘A လြမ္းEတြကUI

၈ဏန္းEတြနဲ႔

သက္Eသၿပရရင္ၿဖင့္ ကUIတUI႔Eနတဲ့ ကမာၻႀကီးဟာ ပန္းကမာၻႀကီး ၿဖစ္သြားပါလIမ့္မယ္ ’ လIU႔ပါ။ ၈ဏန္းEတြနဲ႔ သက္EသၿပရာကEန ဘာလUI႔ ပန္းကမာၻႀကီး ၿဖစ္သြားမလဲဆUIEတာ့ ကၽြန္Eတာ္တUI႔ တြက္Eနတဲ့ ၈ဏန္းEတြဟာ ကမာၻႀကီးသာယာလွပဖUI႔ A သUံးၿပUမွာ ၿဖစ္တဲ့A တြက္ ‘ပန္း’ EလးEတြပါဘဲ လUI႔ Eရးရင္း ဒီ post ကUI A ဆUံးသတ္လUIက္ပါတယ္။

။ Mtssnrty