2014-01-06 3 views
1

MATLAB에서 Newton-Raphson을 반복해야합니다. 그것은 쉬운 것처럼 보이지만 내가 틀린 곳을 알아낼 수는 없습니다. 문제는 : 음 들어값을 변경하고 함수를 수행하는 방법 Newton-Raphson

= 1 : 1) m = 1 개 테이크 C1 = C2 = C1B 및 1-C1 및 u1,2 (I) 및 p1,2 (I)에 대한 루프 않으면

2) m = 2 인 경우, c1 = c1 + dc 및 c2 = 1-c1을 ​​취하여 u1,2 (i) 및 p1,2 (i)에 대해 새로운 c1 및 c2를 갖는 루프를 수행하십시오.

3) m = 3이면 c1 = (c1 * st (1) - (c1-dc) * st (2))/(st (1) -st (2))를 취하여 새로운 c1과 c2에 대한 루프를 수행합니다.

그런 다음 반복 횟수를 늘리십시오. mmm = 2; mmm은 N-R 반복 횟수를 유지합니다. 첫 번째 반복에는 mmm = 1, 두 번째 mmm = 2 등이 있습니다 (이 특정 실행에는 두 번의 반복 만 수행됨). sumint는 내부에 있습니다.

이 수치를 코드에 표시해야하지만 MATLAB은 아래에 오류를 제공합니다. 도와주세요.

코드

관련 부분 : 유 압력의 C1, C2는 캠버 효과 H1DDOT 및 ADDOT 번째 H1의 유도 및 A.의 SUM1되어 있고 SUM2 안쪽되는 속도의 P이다

ii=101; 
    u = cell(2, 1); 
ini_cond = [0,0]; 
for i = 1:2; 
    u{i} = zeros(1,ii); 
    u{i}(:, ii) = ini_cond(i) * rand(1, 1); 
end  


for i=1:ii; 
     fikness=fik*sin(pi.*x); 
     u{1}(i)=(c1-H1D*(x-0.5)+AD/2.*(x-0.5).^2)./(H1-0.5*fikness-A*(x-0.5)); 
     u{2}(i)=(c2+H1D*(x-0.5)-AD/2.*(x-0.5).^2)./(1.-H1+0.5*fikness+A*(x-0.5)); 
end 

p = cell(2, 1); 
q = cell(2, 1); 


for i = 1:2; 
    p{i} = zeros(1,ii); 
    q{i} = zeros(1,ii); 
end  

p{1}(1)=0.5*(1.-u{1}(1).^2); 
q{1}(1)=0; 
p{2}(1)=0.5*(1.-u{2}(1).^2); 
q{2}(1)=0; 

for i=2:101 

q{1}(i)=q{1}(i-1)-dx*(u{1}(i-1)-ub{1}(i-1))./dt; 
p{1}(i)=0.5*(1.-u{1}(i).^2)+q{1}(i); 
q{2}(i)=q{2}(i-1)-dx*(u{2}(i-1)-ub{2}(i-1))./dt; 
p{2}(i)=0.5*(1.-u{2}(i).^2)+q{2}(i); 

end 


st = zeros(2, length(t)); 
st(1,:)=p{1}(100)-p{2}(100); 
m=m+1; 

if m==3; 
c1=(c1*st(1)-(c1-dc)*st(2))/(st(1)-st(2)); 
c2=1-c1; 
end 
for i = 1:2; 
    sumint{i} = zeros(1,length(t)); 
end 
sumint = cell(2, 1); 
sumint{1}(1)=0.5*(p{2}(1)-p{1}(1)); 
sumint{2}(1)=0.5*(p{2}(1)-p{1}(1)).*(-1/2); 

for i=2:ii-1; 
x=(i-1)*dx; 
sumint{1}(i)=sumint{1}(i-1)+(p{2}(i)-p{1}(i)); 
sumint{2}(i)=sumint{2}(i-1)+(p{2}(i)-p{1}(i))*(x-1/2); 

end 

H1DDOT=-sumint{1}.*dx./rmass; 
H1D=H1D+dt*H1DDOT; 
H1=H1+dt*H1D; 
ADDOT=sumint{2}*dx./rmomi; 
AD=AD+dt*ADDOT; 
A=A+dt*AD; 

H1L=H1+A.*0.5; 
H1R=H1-A.*0.5; 
H2=1.-H1; 
rat1=AD./ADinit; 
rat2=ADDOT./AD; 

H1DDOT 및 ADDOT의 값을 정의하기 위해 적분 값을 계산합니다. H1DDOT 및 ADDOT은 시간의 함수입니다. 당신이 메시지에서 볼 수 있듯이

+0

내 목표는 변수 c1, c2를 어떻게 바꾸고 모든 u1,2 (i)와 압력 p1,2 (i) ss에 넣은 다음 H1D, H1DDOT 및 AD, ADDOT .. 고맙습니다. – Meva

답변

0

, 오류가이 라인이다 :

sumint{2}(i) = ... 

이 부분은 당신이 무엇이든 삽입 할 의미 : 이제

sumint{2}(i)=sumint{2}(i-1)+(p{2}(i)-p{1}(i)).*(x-1/2); 

,의는 이유를 알아 보자 셀 sumint{2}에 배열의 오른쪽에있는 i 번째 위치로 이동합니다. 즉, 오른쪽은 스칼라이어야합니다.

?

글쎄, sumint{2}(i-1)+(p{2}(i)-p{1}(i))은 모든 벡터/배열에 대한 색인으로 단일 값을 사용하기 때문에 확실히 스칼라입니다. 문제는 곱셈 .*(x-1/2);입니다.

위의 코드에서 x은 벡터/배열입니다 (length(x) 등을 사용하므로). 스칼라 sumint{2}(i-1)+(p{2}(i)-p{1}(i))에 벡터 x을 곱하면 위에서 설명한대로 작동하지 않는 벡터를 다시 얻을 수 있습니다.

i 값이 x일까요? 코드에서 일어나는 여러 가지 다른 이상한 것들, 예를 들어,이

있습니다

for i=1:101; 
     fikness=fik*sin(pi.*x); 
     u{1}=(c1-H1D*(x-0.5)+AD/2.*(x-0.5).^2)./(H1-0.5*fikness-A*(x-0.5)); 
     u{2}=(c2+H1D*(x-0.5)-AD/2.*(x-0.5).^2)./(1.-H1+0.5*fikness+A*(x-0.5)); 
end 

왜 여기에 루프를해야합니까? 당신은 아무 것도 루핑하지 않고 같은 계산을 100 번하고 있습니다. 다시, 나는 그것이 x(i)이어야한다고 생각한다.

업데이트 : 당신은 스칼라에 x를 변경하기 때문에

새 오류가 도입된다. 그러면 u = zeros(1,length(x))은 스칼라 일 뿐이며 u{1}(i)i ~= 1에 실패합니다.

+0

대단히 감사합니다. 나는 코드를 변경 x는 더 이상 vecto되지 않습니다 :)하지만 matlab에 그것을 실행하면 : ??? u. % cell (2)에 액세스하려고 시도했습니다. numel (u. % 셀) = 1이기 때문에 범위를 벗어납니다. 문제를 해결하는 방법은 무엇입니까? 에있는 ==> grab3_SmithWilson의 오류 p {1} (i) = 0.5 * (1. -u {1} (i).^2) + q {1} (i); – Meva

+0

@ user3115779 : 업데이트를 참조하십시오. 나는 여러분이'x'를 벡터로 유지하기를 정말로 확신합니다. –

+0

네가 맞다.하지만 나는 이미 유속을 u = 0 (1,101)으로 변경하고 모든 루프에 x = (i-1) * dx 라인을 추가했다. 따라서 x는 sclara이고 모든 반복에서 dx의 양과 함께 변경됩니다. 그러나 다시 같은 오류. 나는 분명히 이해하지 못합니다 .. – Meva

관련 문제