2013-04-15 3 views
1

세 가지 미분 방정식 시스템을 모델링하려고합니다. 이것은 방울 모델이며, 호 길이에 대해 매개 변수화 된 것입니다. 입니다.ODE 시스템, 차동 초기 조건이있는 IVP

방정식은 :
DX/DS = COS (세타)
DZ/DS = 죄 (세타)
(세타)/DS = 2 * B + C * Z-죄 (세타)/X

초기 조건은 s = 0에서 x, z 및 theta가 모두 0이라는 것입니다. d (theta)/ds에 특이점을 피하기 위해, 나는 또한 s = 0에 d (theta)/ds = b라는 조건을 가지고 있습니다. 나는 이미 다음과 같은 코드를 작성했다 :

[s,x]=ode23(@(s,x)drpivp(s,x,p),sspan,x0); 

%where p contains two parameters and x0 contains initial angle theta, x, z values. 
%droplet ODE function: 
function drpivp = drpivp(s,x,p); 
%x(1)=theta 
%x(2)=x 
%x(3)=z 
%b is curvature at apex 
%c is capillarity constant 
b=p(1); 
c=p(2); 
drpivp=[2/p(1)+p(2)*x(3)-sin(x(1))/x(2); cos(x(1)); sin(x(1))]; 

어떤 것이 나선형의 해결책을 만들어 낸다. 하나의 작은 물방울 프로파일을 생성하는 대신 많은 물을 생성합니다. 물론, 여기에서는 s = 0에서 theta에 대해 다른 방정식을 사용하는 방법이 확실하지 않기 때문에 방정식을 올바르게 초기화하지 않았습니다.

그래서 질문은 : d (theta)/ds = b가 s = 0 일 때가 아니라 초기 조건을 어떻게 포함시킬 수 있습니까? 이 가능한 matlab에 내장 된 솔버를 사용하고 있습니까?

감사합니다. 이 일을 여러 가지 방법이

답변

1

이며, 가장 쉬운 방법은 단순히 당신의 식으로 if 문을 경우를 추가하는 것입니다

function drpivp = drpivp(s,x,p); 
%x(1)=theta 
%x(2)=x 
%x(3)=z 
%b is curvature at apex 
%c is capillarity constant 
b=p(1); 
c=p(2); 
if (s == 0) 
    drpivp=[b; cos(x(1)); sin(x(1))]; 
else 
    drpivp=[2/p(1)+p(2)*x(3)-sin(x(1))/x(2); cos(x(1)); sin(x(1))]; 
end 
+0

간단하고 효과적인 솔루션! ode 솔버가 그것을 받아 들일 것이라고 생각하지 않았습니다. 감사! (호기심에서 벗어나는 다른 방법을 간략하게 설명해 주시겠습니까?) – Rstallard

+1

@Kud이 문제를 해결하는 또 다른 방법은 시스템을 d2Theta/dt2로 표현하는 것입니다. 이는 차등을 확장해야한다는 것을 의미합니다. – Rasman