Solving a problem with the Shooting Newton Method in Matlab

Question

Francesca 2023-8-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2008202-solving-a-problem-with-the-shooting-newton-method-in-matlab

移动：John D'Errico 2023-9-26

I tried to implement the Shooting method using Newton method to solve a differential equation, but the exact solution recommended by the professor does not coincide with the numerical solution.

Equation:

-y''+exp(-x)/2(y'^2+y^2)=exp(-x)/2*cos(x)+2sin(x)

BC:

2y(0)-y'(0)=0

y(pi)=exp(pi)

Exact solution:

y=exp(x)+sin(x)

The graphs should coincide, but they don't.

Can you help me?

For testing

>> [x,u]= NonLinearShootingNewtonExercise(20,1,0.0001);

>> y=exp(x)+sin(x);

>> plot(x,u,'r*',x,y,'k-')

My code is:

[x,u]= NonLinearShootingNewtonExercise(20,1,0.0001);

y=exp(x)+sin(x);

plot(x,u,'r*',x,y,'k-')

function [x,u] = NonLinearShootingNewtonExercise(N,sk,tol)

%Shooting method using Netwod method

%Equation:

%-y''+exp(-x)/2(y'^2+y^2)=exp(-x)/2*cos(x)+2sin(x)

%BC:

%2y(0)-y'(0)=0

%y(pi)=exp(pi)

%Exact solution:

%y=exp(x)+sin(x)

%Input

%N: number of discretization point

%s0: initial value of s

%tol: tollerance

% Output:

% x: Vector containing the points at which the solution is approximated

% u: Approximated solution

%For testing

% [x,u]= NonLinearShootingNewtonExercise(20,1,0.0001);

% y=exp(x)+sin(x);

% plot(x,u,'r*',x,y,'k-')

epi=exp(pi);

h=pi/N;

x=0:h:pi;

ds=1+tol;

while ds>=tol

[xx,uvUV]=ode23(@fun,x,[sk 2*sk 1 2]); %I have 4 initial values (u,v,U,V)

%uvUV is an array of 4 coloumns and N+1 rows

phik=uvUV(N+1,1)+uvUV(N+1,2)-epi;

dphik=uvUV(N+1,3)+uvUV(N+1,4); %derivate of phi

skp1=sk-phik/dphik;

ds=abs(skp1-sk);

sk=skp1;

end

u=uvUV(:,1);

end

function uvp= fun(x,uvUV)

uvp=[uvUV(2); (exp(-x)/2*(uvUV(2)^2+uvUV(1)^2))-(exp(-x)/2+cos(x)+2*sin(x));...

uvUV(4);(exp(-x)*uvUV(1)*uvUV(3))+(exp(-x)*uvUV(2)*uvUV(4))];

%u is the first component and v is the second component

end

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Francesca 2023-9-26

移动：John D'Errico 2023-9-26

在 MATLAB Online 中打开

%Esercizio n.1 del 22 ottobre

function [x, u] = NonLinearShootingNewtonExercise(N, s0, tol)  
% Shooting method using Newton's method
% Equation:
% -y''+exp(-x)/2(y'^2+y^2)=exp(-x)/2+cos(x)+2sin(x)
% BC:
% 2y(0)-y'(0)=0
% y(pi)=exp(pi)
% Exact solution:
% y=exp(x)+sin(x)
% Input
% N: number of discretization points
% s0: initial value of s
% tol: tolerance
% Output:
% x: Vector containing the points at which the solution is approximated
% u: Approximated solution
% For testing
% [x, u] = NonLinearShootingNewtonExercise(20, 1, 0.0001);
% y = exp(x) + sin(x);
% plot(x, u, 'r*', x, y, 'k-')
esponenzialepi = exp(pi);
h = pi / N;
x = 0:h:pi;
ds = 1 + tol;
sk = s0;
while ds >= tol
    [xx, uvUV] = ode23(@fun, x, [sk 2*sk 1 2]);
    phik = uvUV(N+1, 1)  - esponenzialepi;
    dphik = uvUV(N+1, 3) ; % derivative of phi
    skp1 = sk - phik / dphik;
    ds = abs(skp1 - sk);
    sk = skp1;
end
u = uvUV(:, 1);
u(end) = exp(pi); % Set the boundary condition at x=pi
end
function uvp = fun(x, uvUV)
u = uvUV(1);
v = uvUV(2);
du = uvUV(3);
dv = uvUV(4);
uvp = [v; ...
    (exp(-x) / 2 * (v^2 + u^2)) - exp(-x) / 2 - cos(x) - 2 * sin(x); ...
    dv; ...
    (exp(-x) * u * du) + (exp(-x) * v * dv)];
end
 

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2023-8-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2008202-solving-a-problem-with-the-shooting-newton-method-in-matlab#answer_1287302

编辑：Torsten 2023-8-14

在 MATLAB Online 中打开

sol = fsolve(@fun,23,optimset('TolFun',1e-12,'TolX',1e-12));

Equation solved, inaccuracy possible. The vector of function values is near zero, as measured by the value of the function tolerance. However, the last step was ineffective.

fun(sol)

ans = -1.1977e-12

[X,Y] = ode_solver(sol);

figure(1)

hold on

plot(X,Y(:,1))

plot(X,exp(X)+sin(X))

hold off

grid on

sol = fsolve(@fun,30,optimset('TolFun',1e-12,'TolX',1e-12));

Equation solved, solver stalled. fsolve stopped because the relative size of the current step is less than the value of the step size tolerance squared and the vector of function values is near zero as measured by the value of the function tolerance.

fun(sol)

ans = 3.0909e-13

[X,Y] = ode_solver(sol);

figure(2)

hold on

plot(X,Y(:,1))

plot(X,exp(X)+sin(X))

hold off

grid on

function res = fun(p)

[X,Y] = ode_solver(p);

res = 2*Y(end,1)-Y(end,2);

end

function [X,Y] = ode_solver(p);

fun_ode = @(x,y)[y(2);-exp(-x)/2*cos(x)-2*sin(x)+exp(-x)/2*(y(2)^2+y(1)^2)];

y0 = [exp(pi);p];

tspan = [pi,0];

[X,Y] = ode15s(fun_ode,tspan,y0,odeset('RelTol',1e-12,'AbsTol',1e-12));

end

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Francesca 2023-9-20

How can I modify the code while maintaining the structure and without changing my commands too much?

请先登录，再进行评论。

Solving a problem with the Shooting Newton Method in Matlab

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Solving a problem with the Shooting Newton Method in Matlab

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

WeChat

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论