Energy method with base is based by piecewise linear functions
3 次查看(过去 30 天)
显示 更早的评论
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1530182/image.jpeg)
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 个评论
回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!