Oscillatory solutions problem ODE45
4 ビュー (過去 30 日間)
古いコメントを表示
Hi all,
I'm solving this set of ODEs:
function Y=b(zeta,X)
dF1dzeta=X(2);
dF2dzeta=X(1)^2-X(3)^2+X(2)*X(5)+1;
dG1dzeta=X(4);
dG2dzeta=2*X(1)*X(3)+X(4)*X(5);
dH1dzeta=-2*X(1);
Y = [dF1dzeta; dF2dzeta; dG1dzeta; dG2dzeta; dH1dzeta];
Along with the boundary conditions X(1)=X(3)=X(5)=0 at zeta = 0
and X(1)=0, X(3)=1 at zeta = inf
Normally I would use a shooting method such as this
clear all;
dF2=0.0001;
dG2=0.0001;
b1=30;
initF2=-0.9420;
initG2=0.7729;
K=zeros(2);
zetaspan=linspace(0,b1,b1*100+1);
H=[1;1];
options=odeset('AbsTol',1e-8,'RelTol',1e-6);
while max(abs(H))>1e-10,
[zeta,X]=ode45(@b,zetaspan,[0;initF2+dF2;0;initG2;0],options);
n=size(zeta,1);
X1=[X(n,1);X(n,3)-1];
[zeta,X]=ode45(@b,zetaspan,[0;initF2;0;initG2+dG2;0],options);
n=size(zeta,1);
X2=[X(n,1);X(n,3)-1];
[zeta,X]=ode45(@b,zetaspan,[0;initF2;0;initG2;0],options);
n=size(zeta,1);
X3=[X(n,1);X(n,3)-1];
K(1,1)=(X1(1)-X3(1))/dF2;
K(2,1)=(X1(2)-X3(2))/dF2;
K(1,2)=(X2(1)-X3(1))/dG2;
K(2,2)=(X2(2)-X3(2))/dG2;
H=K\-X3;
initF2=initF2+H(1);
initG2=initG2+H(2);
end
figure;
hold all;
plot(zeta,X(:,1),'LineWidth',2);
plot(zeta,X(:,3),'LineWidth',2);
plot(zeta,X(:,5),'LineWidth',2);
hold off;
xlabel('$\zeta$','interpreter','latex','FontSize',30);
h=legend('$F$','$G$','$H$','Location','NorthEast');
set(h,'interpreter','latex','FontSize',30)
axis([0 30 -1 2]);
grid on
box on
guessing values for F'(0) and G'(0) to solve the problem. However the solutions are oscillatory as zeta tends to infinity and the above method breaks down for zeta>6 (here zeta=b1). So I reformulated the problem, this time shooting from infinity to zero.
clear all;
dF1=0.00001;
dG1=0.00001;
dH1=0.00001;
b1=30;
initF1=0;
initG1=1;
initH1=1.3494;
K=zeros(3);
zetaspan=linspace(b1,0,b1*100+1);
H=[1;1;1];
options=odeset('AbsTol',1e-8,'RelTol',1e-6);
while max(abs(H))>1e-10,
[zeta,X]=ode45(@b,zetaspan,[initF1+dF1;0;initG1;0;initH1],options);
n=size(zeta,1);
X1=[X(n,1);X(n,3);X(n,5)];
[zeta,X]=ode45(@b,zetaspan,[initF1;0;initG1+dG1;0;initH1],options);
n=size(zeta,1);
X2=[X(n,1);X(n,3);X(n,5)];
[zeta,X]=ode45(@b,zetaspan,[initF1;0;initG1;0;initH1+dH1],options);
n=size(zeta,1);
X3=[X(n,1);X(n,3);X(n,5)];
[zeta,X]=ode45(@b,zetaspan,[initF1;0;initG1;0;initH1],options);
n=size(zeta,1);
X4=[X(n,1);X(n,3);X(n,5)];
K(1,1)=(X1(1)-X4(1))/dF1;
K(2,1)=(X1(2)-X4(2))/dF1;
K(3,1)=(X1(3)-X4(3))/dF1;
K(1,2)=(X2(1)-X4(1))/dG1;
K(2,2)=(X2(2)-X4(2))/dG1;
K(3,2)=(X2(3)-X4(3))/dG1;
K(1,3)=(X3(1)-X4(1))/dH1;
K(2,3)=(X3(2)-X4(2))/dH1;
K(3,3)=(X3(3)-X4(3))/dH1;
H=K\-X4;
initF1=initF1+H(1);
initG1=initG1+H(2);
initH1=initH1+H(3);
end
figure;
hold all;
plot(zeta,X(:,1),'LineWidth',2);
plot(zeta,X(:,3),'LineWidth',2);
plot(zeta,X(:,5),'LineWidth',2);
hold off;
xlabel('$\zeta$','interpreter','latex','FontSize',30);
h=legend('$F$','$G$','$H$','Location','NorthEast');
set(h,'interpreter','latex','FontSize',30)
axis([0 30 -1 2]);
grid on
box on
This time using the known values of F(inf)=0 and G(inf)=1 whilst guessing a value for H(inf). But this doesn't seem to get round the problem. What am I best advised to do here? Is there an alternative ODE solver that will combat these oscillatory solutions?
0 件のコメント
回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Ordinary Differential Equations についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!