Energy method with base is based by piecewise linear functions
4 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu Other Formats finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!