MATLAB：調試牛頓拉夫遜法

Question

我試圖建模的液麪控制，我試圖解決的隱方程使用牛頓拉夫遜法

t=(0:1:10); 
x=[14 -4.32E-4 5.28E-4]; y=[14 7]; 
[t,x]=ode45(@PHandLiquidlevelmain,t,x) % returns x 

% For this time period t ,I have x as a 11X3 double array. 
% And I am trying to access that array to calculate y 

y=zeros(max(size(t)),2); 

cv=8.75; 
pk1=6.35; 
pk2=10.25; 


for l=1:1:1% For output Y(1) 
    for k=1:1:11 
     if l<2 
     y(k,l)=y(k,l)+x(k,l);% Equation governed to calculate y(1) 
     end 
    end 
end 

for l=2:1:2 
    for k=1:1:11 

     if l<3 

     syms z; 
     format long; 

     Equa_0=((x(k,l))+(10^(z-14))+((x(k,l+1))*((1+2*(10^(z-pk2)))/(1+(10^(pk1-z))+ 
                (10^(z-pk2)))))-(10^(-z))) 
     % using looping to get values of x(1,2)to x(t,2) and x(1,3) up to x(t,3) 
     % and using it in Equa_0 to create t equations which are need to be 
     % solved using Newton Raphson method 

     %Newton Raphson method  

     e = 1e-5; % setting the tolerance value 
     dx = e + 1; 
     guess = 7; % initially assumed value of z 

     count = 0; % setting counter to know the 
         % no of iterations taken      
     p = zeros(1,1); 
     while (abs(dx) > e) % initialising the iteration and 
          % continue until the error is less than tolerance 

      dx = (eval(Equa_0/(diff(Equa_0)))); % calculating dx, diff is used for 
               % finding the differentiation of the 
               % function 
      guess = guess - dx; % updating the value of x 
      count = count + 1; % incrementing the counter 
      p(count) =guess; 
      drawnow(); 
      plot(abs(p),'r','linewidth',3); 
      grid; 
      if (count > 300) 

       fprintf('Error...! Solution not converging !!! \n'); % printing the 
                     % error message 
       break; 
      end 
     end 
     if (count < 300) 
      fprintf('The solution = '); %printing the result 
      y(k,l)=y(k,l)+guess;% y(1,1) to y(t,1) is calculated 
      fprintf('\nNumber of iteration taken = %d\n',count); 
     end 
     end 
    end 
end 
y 

輸出爲：

t = 

0 
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 


x = 

14.0000 -0.0004 0.0005 
14.0001 -0.0004 0.0005 
14.0001 -0.0004 0.0005 
14.0002 -0.0004 0.0005 
14.0002 -0.0004 0.0005 
14.0003 -0.0004 0.0005 
14.0003 -0.0004 0.0005 
14.0003 -0.0004 0.0005 
14.0004 -0.0004 0.0005 
14.0004 -0.0004 0.0005 
14.0005 -0.0004 0.0005 


Equa_0 = 

10^(z - 14) - 1/10^z + (33*(2*10^(z - 41/4) + 1))/(62500*(10^(z - 41/4) + 10^(127/20 - 
z) + 1)) - 27/62500 

而發生以下錯誤：

The following error occurred converting from sym to double: 
Error using mupadmex 
Error in MuPAD command: DOUBLE cannot convert the input expression into a double array. 

If the input expression contains a symbolic variable, use the VPA function instead. 

Error in PHmain (line 55) 
      p(count) =guess; 

我試過vpa函數但無濟於事。 Equa_0接受x值，但似乎創建了不同的表達式。

Answer 1

一，您指定guess並操縱它，但不要將它推入Equa_0。在計算循環中的dx之前，您需要在每次迭代過程中明確說出z = guess;。

二，一般來說，你不應該使用這個符號工具箱。如果x是恆定的，使用匿名函數（這將需要重寫一些代碼的）：

Equa_0=@(z) ((x(k,l))+(10^(z-14))+((x(k,l+1))*((1+2*(10^(z-pk2)))/(1+(10^(pk1-z))+ (10^(z-pk2)))))-(10^(-z))); 

編輯： 我要指出，你還有其他的事情在你的代碼來解決，但我會離開，部分讓你弄清楚。