Boundary value problem in MATLAB, issue with fsolve

3 次查看（过去 30 天）

carlos g 2020-2-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/505037-boundary-value-problem-in-matlab-issue-with-fsolve

评论： darova 2020-2-12

ws.mat

I am trying to solve and ODE with boundary conditions at `0` and `end`, ws.mat can be downloaded at https://gofile.io/?c=RYzPsO

 clear all
    clc
    global omega eta_max_ode eta F T Tp
    
    eta_max_ode = 10;
    omega=0.76;
    
    load('ws','T','Tp','F','eta')
    IC=[0,1];
    
    opt = optimset('Display','off','TolFun',1E-20);
    FI = fsolve(@(FI) eval_boundary_UVWTI(FI),IC,opt)
    [eta_ode_i, fg_ode_i] = solve_UVWTI(FI);
    
    sol_i = [fg_ode_i];
    eta_i = eta_ode_i;
    plot(eta_i,sol_i(:,1))
    
    function [dfi]=UVWTI(t,fi)
    global omega eta_max_ode eta F T Tp
    
    for i=1:length(eta)
        if eta(i)<t
        else
            inde=i;
            break
        end
    end
    
    A11=0;
    A12=1;
    
    A21=0;
    A22=-(T(inde)^(1-omega))*F(inde)-omega*Tp(inde)/T(inde)+Tp(inde)/T(inde);
    
    dfi=zeros(2,1);
    dfi(1)=A11*fi(1)+A12*fi(2);
    dfi(2)=A21*fi(1)+A22*fi(2);
    
    end
    
    function [eta_ode_i, fg_ode_i] = solve_UVWTI(FI)
    global eta_max_ode eta
    
    options = odeset('RelTol',1e-9,'AbsTol',1e-9);
    
    [eta_ode_i, fg_ode_i] = ode45(@UVWTI,eta,FI,options); 
    size(fg_ode_i);
    end
    
    function [gi] = eval_boundary_UVWTI(FI)
    % Get the solution to the ODE with inital condition F
    [eta_ode_i, fg_ode_i] = solve_UVWTI(FI);
    % Get the function values (for BCs) at the starting/end points
    
    w_start = fg_ode_i(1,1); %w(0) = 0
    w_end = fg_ode_i(end,1); %w(inf) = 0
    
    % Evaluate the boundary function
    gi = [
         w_start
         w_end - 1
         ];
    end

I obtained the correct behaviour this is the solution tending to a constant. However, I should get sol_i(:,1) tending to 1 and fsolve does not seem to calculate the correct initial condition so that this happen. What's wrong in the code? eval_boundary_UVWTI() seems to be correct