Cannot find explicit solution

Question

0 votes

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

13 Comments
Show 11 older comments Hide 11 older comments

Hamza on 3 May 2024

Edited: Torsten on 3 May 2024

Open in MATLAB Online

@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

%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

%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

%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

%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

%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

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 on 3 May 2024

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.

Cannot find explicit solution

13 Comments
Show 11 older comments Hide 11 older comments

Answers (0)

Categories

Products

Release

Tags

Community Treasure Hunt

Cannot find explicit solution

13 Comments Show 11 older comments Hide 11 older comments

Answers (0)

Categories

Products

Release

Tags

See Also

Community Treasure Hunt

13 Comments
Show 11 older comments Hide 11 older comments