Bisection implementation in optimizing non linear ode solution

조회 수: 5 (최근 30일)
Roi Bar-On
Roi Bar-On 2021년 2월 8일
편집: Roi Bar-On 2021년 2월 8일
Hey guys. I have this ODE with the following bc's:
My code includes the shooting method and an intellectual guess of the derivative of rho at x=0 (transforming one BVP problems to one IVP problem), using the Implicit Euler method to write down the derivatives, solve the algebraic system with the Newton-Raphson method and make sure that I converge with infinity norm condition. I try to optimize my results by conducting bisection on the difference between a current value of rho to the its value at the bc that I chose (minimizing f1). This is the implicit Euler code:
function[x,y,f1,BISEC]=implicit_euler_method_DFT(y10,y20,x0,xf,dx,m,n)%i.g=-0.4394
x=x0:dx:xf;
y=zeros(2,length(x));
y(1,1)=y10;
y(2,1)=y20;
B=1;k=1;bc=0.25;
for i=1:1:length(x)-1
y(:,i+1)=y(:,i);
J=Jacovian(y(:,i));
g=[y(1,i+1)-dx*y(2,i+1)-y(1,i);
y(2,i+1)-dx*(2*B/k*y(1,i+1)+1/k*log(abs(y(1,i+1)))-B/k)-y(2,i)];
d=gausselemination(J,g);
f1=y(1,i+1)-bc;
BISEC=bisectionDFTbc(f1);
while (infinity_norm(d)>10^(-m))||(infinity_norm(g)>10^(-n))
y(:,i+1)=y(:,i+1)+d;
J=Jacovian(y(:,i+1));
g=[y(1,i+1)-dx*y(2,i+1)-y(1,i);
y(2,i+1)-dx*(2*B/k*y(1,i+1)+1/k*log(abs(y(1,i+1)))-B/k)-y(2,i)];
d=gausselemination(J,g);
f1=y(1,i+1)-bc;
BISEC=bisectionDFTbc(f1);
end
end
end
where x is the x tilde basically, y is rho tilde, y10 is the initial value for rho, y20 is the initial guess (-0.4394) for its derivative at x=0,x0=0, xf can be 1 or more, dx is to our desire, m and n are tolerances. The 2nd code for the bisection is:
function [x,iter_num]= bisectionDFTbc(a,b,n,m)
iter_num=0;%iteration number
x=(a+b)/2;
while sqrt((a-b)^2)>10^(-m)||sqrt(f1(x)^2)>10^(-n)
iter_num=iter_num+1;
x=(a+b)/2;
if f1(a)*f1(x)<0
b=x;
else
a=x;
end
end
x=(b+a)/2;
end
I have two main problems:
1.The y values go to infinity very fast - I don't get a stable solution.
2.Obviously I don't connect the two functions well.
I would appreciate your help.
Thanks,
Roi

답변 (0개)

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by