Hi ,i have problem whit the function quad
13 views (last 30 days)
Show older comments
% -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)
2 Comments
Geoff Hayes
on 20 Apr 2020
Elkalai - please copy and paste the full error message to this question (rather than creating tags of the error message).
Answers (1)
See Also
Categories
Find more on Object Containers in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!