如何更改值并执行 Newton-Raphson 函数

Question

我需要在 MATLAB 中迭代 Newton-Raphson。这似乎很容易，但我无法弄清楚我错在哪里。问题是：

对于 mmm=1： 1) 如果 m=1 取 c1=c1b 和 c2=1-c1 并为 u1,2(i) 和 p1,2(i) 执行循环

2)如果 m=2 取 c1=c1+dc 和 c2=1-c1，这次用新的 c1 和 c2 为 u1,2(i) 和 p1,2(i) 循环

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，依此类推（此特定运行仅进行 2 次迭代）。 sumint 在积分内。

我需要在代码中绘制这些数字，但 MATLAB 在下面给出错误。请帮帮我。

相关部分代码：

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;

u 是速度 p 是压力 c1,c2 是外倾效应 H1DDOT 和 ADDOT 是 H1 和 A 的二阶导数。 sum1 和 sum2 是定义 H1DDOT 和 ADDOT 值的积分内部。 H1DDOT 和 ADDOT 是时间的函数。

Answer 1

从消息中可以看出，错误在于这一行：

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

现在，让我们找出原因：

sumint{2}(i) = ...

这部分意味着您想要将右侧的任何内容插入到单元格sumint{2} 中数组的ith 位置。这意味着，右侧必须是标量。

是吗？

好吧，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，将返回一个向量，如上所述，这不起作用。

也许您想要x 的ith 值？

您的代码中还发生了其他一些奇怪的事情，例如：

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 将失败。