Energy method with base is based by piecewise linear functions

3 次查看（过去 30 天）

显示更早的评论

Aurora 2023-11-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2043177-energy-method-with-base-is-based-by-piecewise-linear-functions

评论： Aurora 2023-11-5

在 MATLAB Online 中打开

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