Cannot find explicit solution

Hiii ,i need some help in matlab concerning programming pseudo Spectral method for Volterra equation i have all the script but when i check the solution i found i message error Warning: Cannot find explicit solution. > In solve (line 316) In HMM (line 350) La valeur numérique de x lorsque t = 0.0625 est :

can you help me please

[Torsten]

Use a numerical, not a symbolic approach.

[Hamza]

how I use it, I'm new on matlab

[Sam Chak]

Which other computational software have you utilized for solving numerical problems? The logical skills required for problem-solving are generally transferrable to MATLAB environments when addressing numerical problems.

[Hamza]

I am using matlab2016a, can I send you my script and you can solve this issue for me?

I really need this programme, I have a graduation next month and I need it in my thesis.

[Torsten]

If you are new to MATLAB, it's too late if you need it next month.

Maybe you can find MATLAB code in the internet where the Pseudospectral method for differential equations is already implemented.

[Walter Roberson]

It is possible that there is no symbolic solution.

[Sam Chak]

@Hamza, maintain your determination; your goal is within reach already. Encourage yourself to accomplish the task in the forthcoming days. Assuming you possess a mathematical background (which I presume you do), dare to make conjectures enabling the transformation of the Volterra equation into a special class of nonlinear equations for which explicit solutions are known.

Subsequently, revise your query accordingly and share the MATLAB script in this forum.

[John D'Errico]

We cannot help you here, if you will not post your code. Answers is not that type of site.

[Hamza]

@John D'Errico

% input variable
b = 1; 
T = b; 
N = 4; 
M = N; 


% D?finir les conditions initiales et aux bords 
syms x t ; 
h0 = exp(x); 
h1 = exp(t); 
h2 = exp(t); 

% D?finir la fonction u(x,t) 
u = @(x,t) exp(x+t); 
%disp('La solution g?n?rale est :'); 
%disp(u(x,t)); 

%d?finir la fonction arbitraire 
u_0 = @(x,t) x.*(1-t).*exp(-x.*t); 
%disp('La fonction arbitraire est :'); 
%disp(u_0(x,t)); 

% D?finir la fonction g(x, t) 
syms g(x, t); 
g(x, t) = exp(2*(x+t)); 

%D?finir les vecteur x_k et x_i
x_k = zeros(N+1, 1);   
for i = 0:N 
    x_k(i+1) = ((1-cos(((2*i-1)*pi))/(2*(N+1))))/2;  
end 
x_i = x_k; 
%disp('Les x_i sont :'); 
%disp(vpa(x_i, 4)); 
x_j = zeros(N+1, 1);  
for j = 0:M 
    x_j(j+1) = ((1-cos(((2*j-1)*pi))/(2*(M+1))))/2;  
end 
%disp('Les x_j sont :'); 
%disp(vpa(x_j, 8)); 

% Calculer si et tj
s_i = (b + b * x_j) / 2; 
%disp('Les valeurs de s_i sont :'); 
%disp(vpa(s_i, 4)); 
t_j = (T + T * x_j) / 2; 
%disp('Les valeurs de t_j sont :'); 
%disp(vpa(t_j, 8));

% Initialiser le vecteur pour stocker les valeurs de g(x, t_j) 
g_tj_values = sym(zeros(N+1, 1)); 

% Boucle sur les ?l?ments de t_j pour calculer g(x, t_j) pour chaque valeur de t_j 
for i = 1:numel(t_j) 
    
    % Substituer la valeur t_j(i) ? la place de t dans g(x, t) 
    g_tj_values(i) = subs(g, t, t_j(i)); 
end 

% Afficher le r?sultat 
%disp('g(x, t_j) :'); 
%disp(vpa(g_tj_values,8));

% Initialisation des vecteurs des conditions aux bords 
h0_vector = zeros(N+1, 1); 
h1_vector = zeros(M, 1); 
h2_vector = zeros(M, 1); 
h0_si=subs(h0,x,x_j); 
h1_tj=subs(h1,t,x_j); 
h2_tj=subs(h2,t,x_j); 
%disp(vpa(h0_si,8)); 
%disp(vpa(h1_tj,8)); 
%disp(vpa(h2_tj,8)); 

% D?finir les polyn?mes de Bernoulli P_{i,b}(x) et P_{j,T}(t) 
P_x = sym(zeros(N+1,1)); % polyn?mes par rapport ? x 
P_t = sym(zeros(M+1,1)); % polyn?mes par rapport ? t 
for i = 0:N 
    sum_term_x = 0; 
    for k = 0:i 
        sum_term_x = sum_term_x + (-1)^(i-k) * nchoosek(i, i-k) * nchoosek(i+k, k) * (x/b)^k; 
    end 
    P_x(i+1) = sqrt((2*i+1)/b) * sum_term_x; 
end 
for j = 0:M 
    sum_term_t = 0; 
    for l = 0:j 
        sum_term_t = sum_term_t + (-1)^(j-l) * nchoosek(j, j-l) * nchoosek(j+l, l) * (t/T)^l; 
    end 
    P_t(j+1) = sqrt((2*j+1)/T) * sum_term_t; 
end 
%disp('Polynomes par rapport a x :'); 
%disp(vpa(P_x, 8));  
%disp('Polynomes par rapport a t :'); 
%disp(vpa(P_t, 8)); 


% Cr?er la matrice T_{N,b} 
T_b = zeros(N+1,N+1); 
% Remplir la partie inf?rieure de la matrice Tb 
for i = 1:N+1 
    for j = 1:i 
        T_b(i,j)= (-1)^(i-j) * sqrt(2*i-1) * (1/b^(j-0.5)) * nchoosek(i-1,i-j) * nchoosek(i+j-2,j-1);
    end 
end
%disp('La matrice T_b :'); 
%disp(vpa(T_b,8)); 

% Cr?er la matrice T_{M,T} 
T_b = zeros(M+1,M+1); 
% Remplir la partie inf?rieure de la matrice Tb 
for i = 1:M+1 
    for j = 1:i 
        T_b(i,j)= (-1)^(i-j) * sqrt(2*i-1) * (1/T^(j-0.5)) * nchoosek(i-1,i-j) * nchoosek(i+j-2,j-1);
    end 
end
%disp('La matrice T_b :'); 
%disp(vpa(T_b,8)); 

% Initialisation du vecteur Psi 
Psix = sym(zeros(N+1, 1)); % Initialisation d'un vecteur de taille N+1 avec des valeurs symboliques 

% Remplissage du vecteur Psi avec les puissances successives de x 
for i = 1:N+1 
    Psix(i) = x^(i-1); % ?lever x ? la puissance correspondante 
end 
%disp('Psi de x est :'); 
%disp(vpa(Psix,8)); 

% Initialisation du vecteur Psi de t 
Psit = sym(zeros(M+1, 1)); % Initialisation d'un vecteur de taille N+1 avec des valeurs symboliques 

% Remplissage du vecteur Psi avec les puissances successives de t 
for i = 1:M+1 
    Psit(i) = t^(i-1); % ?lever x ? la puissance correspondante 
end 
%disp('Psi de t est :'); 
%disp(vpa(Psit,8)); 

%D?finir les vecteur Phi par rapport a x et t et b et T
Phix= T_b * Psix; 
Phit= T_b * Psit; 
%disp('le vecteur Phi de x est :'); 
%disp(vpa(Phix,8)); 
%disp('le vecteur Phi de t est :'); 
%disp(vpa(Phit,8)); 
Phi0 =  subs(Phix,x, 0);
Phib = subs(Phix,x, b);
PhiT = subs(Phit,t, T);
%disp('le vecteur Phi de 0 est :'); 
%disp(vpa(Phi0,8)); 
%disp('le vecteur Phi de b est :'); 
%disp(vpa(Phib,8)); 
%disp('le vecteur Phi de T est :'); 
%disp(vpa(PhiT,8)); 

%D?finir les vecteur Phi par rapport a si et tj
Phi_si=subs(Phix,x,s_i(i));
%disp('le vecteur Phi_si de x est :'); 
%disp(vpa(Phi_si,8)); 
Phi_tj=subs(Phix,x,t_j(j));
%disp('le vecteur Phi_tj de t est :'); 
%disp(vpa(Phi_tj,8)); 

U = sym(zeros(N+1, M+1));  
for i = 1:N+1 
    for j = 1:M+1  
        %approximer l'int?grale en utilisant la m?thode des quadratures 
        integrand = u_0(x,t) * P_x(i) * P_t(j); 
        U(i,j) = int(integrand, x, 0, b) * int(integrand, t, 0, T); 
    end 
end 
%disp('La matrice U est :'); 
%disp(vpa(U,8)); 
la_taille_de_matrice_U = size(U) ;


U_N_M= Phix' * U * Phit; 
%disp('la fonction u_{N,M} APPROROXIMMER est :'); 
%disp(vpa(U_N_M,8)); 

% Cr?ation de la matrice A_N 
A_N = zeros(N+1); 
for i = 1:N+1 
    for j = 1:N+1 
        if i == j+1 
            A_N(i, j) = i - 1; 
        end 
    end 
end 
%disp('La matrice A_N est :'); 
%disp(A_N); 


% Cr?ation de la matrice A_N 
A_M = zeros(N+1); 
for i = 1:M+1 
    for j = 1:M+1 
        if i == j+1 
            A_M(i, j) = i - 1; 
        end 
    end 
end 
%disp('La matrice A_M est :'); 
%disp(A_M);

% Cr?ation de la matrice B_N 
B_N = zeros(N+1); 
for i = 3:N+1 
    B_N(i, i-2) = (i-1)*(i-2); 
end 
%disp('La matrice B_N est :'); 
%disp(B_N); 

% Cr?ation de la matrice B_N 
B_M = zeros(N+1); 
for i = 3:N+1 
    B_M(i, i-2) = (i-1)*(i-2); 
end 
%disp('La matrice B_M est :'); 
%disp(B_M); 

% Calcule la deriver premier de la fonction U_N_M par rapport a t 
dtu = Phix' * U * T_b * A_N * inv(T_b) * Phit;
% Calcule la deriver secend de la fonction U_N_M par rapport a x 
dxxu = (T_b * B_M * inv(T_b) * Phix)' * U * Phit;

%disp('la deriver premier de la fonction U_N_M par rapport a t :'); 
%disp(vpa(dtu,8)); 
%disp('la deriver secend de la fonction U_N_M par rapport a x :'); 
%disp(vpa(dxxu,8)); 

syms x;
% Boucle pour g?n?rer et afficher les polyn?mes de Legendre
for n = 0:N
    legendre_poly = legendreP(n, x);
   % disp(['Polyn?me de Legendre P', num2str(n), '(x) : ']);
   % disp(vpa(legendre_poly,8));
end

% Boucle pour calculer et afficher les racines de chaque polyn?me de Legendre
for n = 0:N+1
    legendre_poly = legendreP(n, x);
    roots = vpasolve(legendre_poly == 0, x);
   % disp(['Les racines du polyn?me de Legendre P', num2str(n), '(x) sont :']);
   % disp(vpa(roots,8));
end

weights = zeros(1, N+1); % Initialiser un vecteur pour stocker les poids

for k = 0:N
    xk = vpasolve(legendreP(N+1, x) == 0, x); % Trouver les racines du polyn?me de Legendre N+1
    Lr2p1_prime = diff(legendreP(N+1, x)); % D?riv?e du polyn?me de Legendre N+1
    
    % ?valuer la d?riv?e et le poids pour la racine k
    Lr2p1_prime_k = subs(Lr2p1_prime, x, xk(k+1));
    weight_k = T / ((1 - xk(k+1))^2 * Lr2p1_prime_k^2);
    
    weights(k+1) = weight_k; % Stocker le poids dans le vecteur des poids
end

%disp('Les poids associ?s aux racines des polyn?mes de Legendre sont :');
%disp(weights);

syms s t; 
y=exp(s-t); 

% Initialisation du vecteur des r?sultats 
yt_values = sym(zeros(size(s_i))); 

% Boucle sur les ?l?ments de s_i et t_j pour ?valuer y pour chaque paire de valeurs 
for i = 1:numel(s_i) 
    % Substitution des valeurs dans y 
    yt_values(i) = subs(y, [s, t], [s_i(i), t_j(i)]); 
end 

% Affichage des r?sultats 
%disp('Les valeurs de yt sont :'); 
%disp(vpa(yt_values, 8)); 

u_0_x_si = u_0(x,s_i); 
%disp('La fonction u_0(x,s_i) est :'); 
%disp(vpa(u_0_x_si, 8)); 


k_1=yt_values' * u_0_x_si; 
%disp('la fonction K_1(x,t_j,s_i,u(x,s_i)) est'); 
%disp(k_1); 

% Ajouter le produit de la fonction k_1 et les poids de Gauss-Legendre ? la somme de l'int?grale 
integral_value = weights * k_1; 
%disp('formule de gauss legendre quadrature est :'); 
%disp(vpa(integral_value',8)); 

% Calcule de la fonction RES(x,t) 
RES= Phix' * U * T_b * A_N * inv(T_b) * Phit  + (T_b * B_N * inv(T_b) * Phix)' * Phit - integral_value' - g_tj_values; 
%disp(vpa(RES,4)); 
taille_de_res=size(RES) ;


part1= Phi_si' * U * T_b * A_N * inv(T_b) * Phi_tj ;
%disp('part1 est'); 
%disp(part1); 
pp=size(part1); 
part2= (T_b * B_N * inv(T_b) * Phi_si)' * U * Phi_tj ; 
%disp('part2 est:'); 
%disp(part2); 
Res_value= part1  + part2 - integral_value' - g_tj_values;
%disp('Res_ value est ');
%disp(vpa(Res_value,8));
taille_de_res=size(Res_value);

eqn1=  Res_value(2:N,:); 
l=size(eqn1) ;

eqn2=Phi_si' * U * Phi0 - h0_si; 
eqn3=Phi0' * U * Phi_tj - h1_tj; 
eqn4= Phib' * U * Phi_tj - h2_tj; 
eqn2=eqn2(1:N+1,:); 
eqn3=eqn3(2:M+1,:);
eqn4=eqn4(2:M+1,:);

%disp('eqn1 est '); 
%disp(vpa(eqn1,8)); 
%disp('eqn2 est '); 
%disp(vpa(eqn2,8)); 
%disp('eqn3 est '); 
%disp(vpa(eqn3,8)); 
%disp('eqn4 est '); 
%disp(vpa(eqn4,8)); 

P=[eqn1;eqn2;eqn3;eqn4];
size(P)

ans = 1x2

16 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

eqn1==0;
eqn2==0;
eqn3==0;
eqn4==0;

t_0 = 0.0625;

% Substituer la valeur de t dans les ?quations
eqn1_substituted = subs(eqn1, t, t_0)

eqn1_substituted = 

eqn2_substituted = subs(eqn2, t, t_0)

eqn2_substituted = 

eqn3_substituted = subs(eqn4, t, t_0)

eqn3_substituted = 

eqn4_substituted = subs(eqn4, t, t_0)

eqn4_substituted = 

% R?soudre pour x
%x_value = solve([eqn1_substituted, eqn2_substituted,eqn3_substituted,eqn4_substituted], x);

% Convertir la solution symbolique en valeur num?rique
%x_numeric = double(x_value); % Convertir en double pour obtenir une valeur num?rique
x_numeric = vpasolve(eqn1_substituted(1)==0,x)

x_numeric = 

% Afficher la valeur de x lorsque t = t_0
disp(['La valeur numerique de x lorsque t = ', num2str(t_0), ' est :']);

La valeur numerique de x lorsque t = 0.0625 est :

disp(x_numeric); 

[Torsten]

eqn1_substituted == 0 , e.g., is a vector of 3 equations in a scalar variable x.

Do you want to solve the three equations separately and get three different x-values ? I did this in your code for eqn1_substituted(1) == 0, thus for the first equation. Maybe there are other solutions apart from x = -104.60...

[Hamza]

@Torsten i want to solve this system the first is eqn1, 2nd is eqn2 ....

to find the approximate solution and at the end i want to calculate the error for differant value of t like this table

Not the values in the table because I used other exemple, but just to give you an idea

[Torsten]

I don't know anything about the mathematical background of your code and won't dive into it - I just wanted to show you that the "solve" command will not work because you have 16 equations in 1 unknown.

So what we can try to do is to find technical errors in your code, but we won't be able to rectify the mathematics.

[Hamza]

@Torsten Thank you anyway