Energy method with base is based by piecewise linear functions
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
The differential problem is:
-(x^2 y')'+6y=8x+1/x
y(-2)-2y'(-2)=-4
y(-1/2)=-5/4
The output isn't little different in te point zero but i don't see the motivation. Can you help me?
function [x,u] = ex_PWLener2(m,n)
%metodo energia con base delle funzioni lineari a tratti con condizione
%essenziale in b
B=11/14;
C=-6/7;
h=1.5/m;
xC=-2:h:-0.5; %x di reazione della base
A=zeros(m);
b=zeros(m,1);
A(1,1)=2+integral(@(x)fA(x,0,0,xC),xC(1),xC(2));
A(1,2)=integral(@(x)fA(x,0,1,xC),xC(1),xC(2));
b(1)=integral(@(x)fb(x,0,xC),xC(1),xC(2));
for k=1:m-2
A(k+1,k+1)=integral(@(x)fA(x,k,k,xC),xC(k),xC(k+2));
A(k+1,k+2)=integral(@(x)fA(x,k,k+1,xC),xC(k+1),xC(k+2));
A(k+1,k)=integral(@(x)fA(x,k,k-1,xC),xC(k),xC(k+1));
b(k+1)=integral(@(x)fb(x,k,xC),xC(k),xC(k+2));
end
A(m,m)=integral(@(x)fA(x,m-1,m-1,xC),xC(m-1),xC(m+1));
A(m,m-1)=integral(@(x)fA(x,m-1,m-2,xC),xC(m-1),xC(m));
b(m)=integral(@(x)fb(x,m-1,xC),xC(m-1),xC(m+1));
delta=A\b;
h=1.5/n;
x=-2:h:-0.5;
u=zeros(1,n+1);
for k=1:m-1
u=u+delta(k)*Cfun(x,k,xC);
end
u=u+B.*x+C;% u=u+l(x) dove l(x) va calcolata
end
function y = Cfun(x,k,xC)
m=length (xC)-1;
y=zeros(size(x));
for j=1:length (x)
if k==0
if xC(1)<=x(j) && x(j)<xC(2)
y(j)=(xC(2)-x(j))/(xC(2)-xC(1));
end
elseif k==m
if xC(m)<x(j) && x(j)<=xC(m+1)
y(j)=(x(j)-xC(m))/(xC(m+1)-xC(m));
end
else

if xC(k)<x(j) && x(j)<=xC(k+1)
y(j)=(x(j)-xC(k))/(xC(k+1)-xC(k));
elseif xC(k+1)<=x(j) && x(j)<xC(k+2)
y(j)=(xC(k+2)-x(j))/(xC(k+2)-xC(k+1));
end
end
end
end
function y=Cder(x,k,xC)
%mettere derivata rispetto x(j)
m=length (xC)-1;
y=zeros(size(x));
for j=1:length (x)
if k==0
if xC(1)<=x(j) && x(j)<xC(2)
y(j)=-1/(xC(2)-xC(1));
end
elseif k==m
if xC(m)<x(j) && x(j)<=xC(m+1)
y(j)=1/(xC(m+1)-xC(m));
end
else
if xC(k)<x(j) && x(j)<=xC(k+1)
y(j)=1/(xC(k+1)-xC(k));
elseif xC(k+1)<=x(j) && x(j)<xC(k+2)
y(j)=-1/(xC(k+2)-xC(k+1));
end
end
end
end
function y=fA(x,k,l,xC)
y=(x.^2).*Cder(x,k,xC).*Cder(x,l,xC)+6.*Cfun(x,k,xC).*Cfun(x,l,xC);
end
function y=fb(x,k,xC)
B=11/14;
C=-6/7;
y=Cfun(x,k,xC).*(8.*x+(1./x)-4*B*x-6*C);
end
2 commentaires
Réponses (0)
Voir également
Catégories
En savoir plus sur Other Formats dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!