Hi ,i have problem whit the function quad

Question

Elkalai Mouncif 2020-4-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/519372-hi-i-have-problem-whit-the-function-quad

回答： darova 2020-4-21

% -u''=(3/2)*t^(1/2)+(2/3)*t^(3/2)+(2/3)-(2/3)*t, int(u'v')=int(fv)
% u(0)=u(1)=0
clear;figure;
n=10;% Nombre des noeuds
a=0;b=1;
x=linspace(a,b,n);h=(b-a)/(n-1);% Noeuds
syms t ;
f=-(3/2)*t^(1/2)+(2/3)*t^(3/2)+(2/3)-(2/3)*t;  
for i=1:n-2
    H1=@(t)(x(i+2)-t)/h;
    H2=@(t)(t-x(i))/h;
    k(i,i)=2/h;
    f1=@(t)f.*(x(i+2)-t)/h;
    f2=@(t)f.*(t-x(i))/h;
    F(i)=quad(f2,x(i),x(i+1))+quad(f1,x(i+1),x(i+2));
end
% Calcul des éléments non diagonaux de la matrice de régidité
for i=2:n-2
    k(i-1,i)=-1/h;
    k(i,i-1)=k(i-1,i);
end
% Matrice globale de régidité
A=zeros(n-2);
for i=1:n-2
    A(i,i)=k(i,i);
end
for i=2:n-2
    A(i-1,i)=k(i-1,i);
    A(i,i-1)=A(i-1,i);
end
%Vecteur de second memebre du système global
B(:,1)=F(1:n-2);
Uh=zeros(n,1); 
Uh(2:n-1)=A\B;% Solution approchée
t=a:.001:b;
u=@(t)2*t.^(3/2)/3-2*t/3 ;% Solution exacte
y=feval(u,t);
plot(t,y),hold on
%m=length(x);
u=@(x)2*x.^(3/2)/3-2*x/3;
U(1:n,1)=feval(u,x(1:n));
Uh(1)=0;Uh(n)=0;
for i=1:n-1
    t=x(i):.001:x(i+1);
    H2=(t-x(i))/h;
    H1=(x(i+1)-t)/h;
    uh=Uh(i)*H1+Uh(i+1)*H2;%Expression analytique de la solution approchée
    plot(t,uh,'-r')
end
%Affichage des résultats
er=abs(Uh-U);
T(:,1)=x';T(:,2)=U;T(:,3)=Uh;T(:,4)=er;
disp(['      x  ','       u ','       uh ','      Erreur'])
disp(T)