猪冰龙

导航

学习newton raphson and back eluer

 1 % % time step  https://ww2.mathworks.cn/matlabcentral/answers/184200-newton-raphson-loop-for-backward-euler
 2 % h = (t_final - t_init)/n; % with n number of time steps
 3 % % vectors
 4 % t = [tinit zeros(1,n)]; % time
 5 % y = [yinit zeros(1,n)]; % solution
 6 % % Backward Euler loop
 7 % for i = 1:n
 8 %    t(i+1) = t(i) + h;
 9 %    y_temp = y(i) + h(f(t(i), y(i)));
10 %    y(i+1) = y(i) + h*f(t(i+1), y_temp);
11 % end
12 % for i = 1:n
13 %    error = 1;
14 %    tolerance = 1e-6;
15 %    t(i+1) = t(i) + h;
16 %    y_temp = y(i) + h*(f(t(i), y(i)));
17 %    while error >= tolerance 
18 %       y(i+1) = y(i) + h*f(t(i+1), y_temp);
19 %       error = abs(y(i+1) - y_temp) % (local) absolute error
20 %       y_temp = y(i+1);
21 %    end
22 % end
23 
24 % yold = y(i)+h*f(t(i),y(i));
25 % while error >= tolerance
26 % ynew = yold-(yold-(y(i)+h*f(t(i+1),yold)))/(1-h*df(t(i+1),yold));
27 % error = abs(ynew-yold);
28 % yold=ynew;
29 % end
30 % y(i+1) = ynew;
31 
32 %y'=y+2*x/y^2 x=[0,2] y(0)=1  https://wenku.baidu.com/view/d18cdaa10b4c2e3f5627632f.html
33 t_final=2;
34 t_init=0;
35 n=5;
36 tolerance=0.0000001
37  h = (t_final - t_init)/n;
38  ti=t_init+h;
39 yold=1+h*f(0,1);% yold = y(i)+h*f(t(i),y(i));
40 while error >= tolerance
41 ynew = yold-(yold-(y(i)+h*f(t(i+1),yold)))/(1-h*df(t(i+1),yold));
42 error = abs(ynew-yold);
43 yold=ynew;
44 end
45 y(i+1) = ynew;

 上面代码应该怎样修改?

posted on 2018-04-25 18:04  猪冰龙  阅读(215)  评论(0编辑  收藏  举报