Hi ,i have problem whit the function quad
13 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
% -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 comentarios
Geoff Hayes
el 20 de Abr. de 2020
Elkalai - please copy and paste the full error message to this question (rather than creating tags of the error message).
Respuestas (1)
Ver también
Categorías
Más información sobre Object Containers en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
