Cannot find explicit solution

Question

Hamza 2024-5-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2114726-cannot-find-explicit-solution

评论： Hamza 2024-5-3

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 个评论
显示 11更早的评论隐藏 11更早的评论

Hamza 2024-5-3

编辑：Torsten 2024-5-3

在 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 2024-5-3

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 2024-5-3