Metodo di collocazione non lineare per risolvere problemi differenziali

il problema è il seguente

-y''+e^-x yy'=e^x (x^2-2)

y(-1)=-1/e

y(1)=e

sol y=xe^x

Non mi sembra di aver sbagliato niente eppure quando run la funzione con m=10, n=2, d0=-0.5 e tol=0.1 non riesce a darmi una risposta perchè gli serve infinito tempo epr calcolare J(l,k+1). Come posso fare per risolvere?

function [x,u] = ex_collNL(m,n,d0,toll)
%due condizioni essenziali
delta = d0*ones(m+1,1); %LASCIA
h = 2/(m+2);  %h=b-a/(m+2)
xC = -1+h:h:3*1-h; %Nodi di collocazione
normphi = toll + 1; %LASCIA
while normphi >= toll   %Inizia il Newton
    J = zeros(m+1);  
    phi = zeros(m+1,1); 
    for l = 1:m+1       %Da qui costruiamo la J, che sta sui fogli.
        for k = 0:m
            %J(elle,k+1) = derivata di epsilon(x_l) rispetto a delta_k.
            J(l,k+1) = -DDu(k,xC(l)) +exp(-xC(l))*(ub(k,xC(l))*Dufun(delta,xC(l))+Du(k,xC(l))*ufun(delta,xC(l)));
        end
        phi(l) = DDufun(delta,xC(l))+exp(-xC(l))*ufun(delta,xC(l))*Dufun(delta,xC(l))-exp(xC(l))*(xC(l)^2-2);
    end
    deltanew=J\phi;
    normphi = sqrt(sum((delta-deltanew).^2));
end
h = 2/n;
x = -1:h:1;
u = ufun(delta,x);
end

%ub
function y = ub(k,x)  %sono le u_k calcolate nel foglio.
y = x.^(k+2)-x.^(k);
end

function y = Du(k,x)  %derivate delle ub precedenti
y = (k+2)*x.^(k+1)-k*(x.^(k-1));
end

function y = DDu(k,x)  %derivate di Di precedenti
y = (k+1)*(k+2)*x.^(k)-(k-1)*k*(x.^(k-2));
end


%ufun
function y = ufun(delta,x)  % ufun = u sui fogli
m = length(delta)-2;
B = (exp(1)^2-1)/(2*exp(1)); %calcolata nei fogli.
A=(exp(1)^2+1)/(2*exp(1));
y = A.*x+B;
for i = 0:m
    y = y + delta(i+1)*ub(i,x);
end
end

function y = Dufun(delta,x) %E' la derivata della u sopra, calcolata sui fogli
m = length(delta)-2; %Stessa della u sopra
A=(exp(1)^2+1)/(2*exp(1));
y = A.*zeros(size(x));
for i = 0:m
    y = y + delta(i+1)*Du(i,x);
end
end

function y = DDufun(delta,x)  %E' la derivata di Dufun = u', calcolata sui fogli
m = length(delta)-2;
y = zeros(size(x));
for i = 0:m
    y = y + delta(i+1)*DDu(i,x);
end
end

Answer 1

[Zuber Khan]

Hi,

I understand that you are trying to solve a boundary value problem (BVP) for a second-order ordinary differential equation (ODE) using a numerical method implemented in MATLAB, and you are encountering issues with the computation taking an infinite amount of time, particularly with the calculation of the Jacobian matrix J(l,k+1).

From the code, it seems you are using some form of collocation method with an iterative Newton-Raphson approach to solve the nonlinear system of equations. The code constructs a Jacobian matrix "J" and a residual vector "phi" to find the corrections "delta" to the approximate solution "u".

As per my understanding, the variable "deltanew" contains the next value of "delta" and you want to upgrade it with every iteration of "while" loop in the "ex_collNL" function. The value of "deltanew" in turn depends on the Jacobian matrix "J" and vector "phi". However, note that the value of "delta" is nowhere updated in the while loop. It is just initialized once in "ex_collNL" function and the same "delta" is used in every iteration while updating the matrix "J" and vector "phi". Kindly observe the snapshot of the respective code lines as follows for a better understanding :

Now, because of the above mentioned reason, the values of "J" and "phi" are not updated correctly and hence the variable "deltanew" is also not converging, leading to an infinite "while" loop. I will recommend to make the necessary changes in the code such that the value of variable "delta" inside the loop is also updated.

One way to do that is to store the initial value of the variable "delta" in another variable, say "deltaold" and then update the value of delta with a new value in every iteration as shown below:

I am getting the following output after making the above changes and the code runs within a few seconds.

I hope it resolves your query.

Best Regards

Zuber