Energy method with Hermite cubic function (2 essential condition)

Question

Aurora 2023-11-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2043132-energy-method-with-hermite-cubic-function-2-essential-condition

评论： Aurora 2023-11-5

采纳的回答： Torsten

在 MATLAB Online 中打开

The different proble is:

-(x^2 y')'+20y=7*x^2 -5

y(-1)=0

y(1)=0

sol: y=(x^2)/2-(x^4)/4-1/4

The next code has the following problem:

Unrecognized function or variable 'm'.

Error in ex_enerHermite>support (line 170)

if l>=2*m

Error in ex_enerHermite (line 10)

[lbl,rbl]=support(l,m);

how can I solve it? Where is the errore?

function [x,u]=ex_enerHermite(m,nv)
%problema con due condizioni essenziali
h=2/m;
xc=-1:h:1;
A=zeros(2*m);%se non ho condizioni essenziali diventa 2*m+2
b=zeros(2*m,1);%se non ho condizioni essenziali diventa 2*m+2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Caso condizione essenziale in a
for l=1:2*m %se non ho condizioni essenziali diventa 2*m+2
    [lbl,rbl]=support(l,m);
    for j=1:2*m %se non ho condizioni essenziali diventa 2*m+2
        if abs(l-j)<4 %altrimenti l'integrale fa zero causa supporto base
          [lbj,rbj]=support(j,xc);
            lb=max(lbl,lbj);
            rb=min(rbl,rbj);
            tol=10^-8;
            if abs(lb-rb)>=tol
% termine di bordo già inserito alla linea 26
                A(l,j)=integral(@(x)fA(x,l,j,xc),lb,rb);
            end
        end
    end
    b(l)=integral(@(x)fb(x,l,xc),lbl,rbl);
end
delta=A\b;
h=2/nv;
x=-1:h:1;
u=zeros(1,size(x));
for l=1:2*m %da cambiare in base a dove fai la sommatoria della soluzione u
    u=u+delta(l)*fu(x,l,xc);
    %se cond. ess. in a allora l=1:2*m+1 e 
    %u=u+delta(l)*fu(x,l,xc)  
    %se non ho cond. ess. allora l=0:2*m+1 e 
    %u=u+delta(l+1)*fu(x,l,xc) 
end
%u Ã¨ la soluzione del problema traslato quindi se necessario aggiungo l(x)
end
function y=fH(x,l,xc) %non toccare
m=length(xc)-1;
y=zeros(size(x));
for j=1:length(x)
    if l==0
        if xc(1)<=x(j) && x(j)<xc(2)
            y(j)=3*(((xc(2)-x(j))/(xc(2)-xc(1)))^2)-...
                2*(((xc(2)-x(j))/(xc(2)-xc(1)))^3);
        end
    elseif l==m
        if xc(m)<x(j) && x(j)<=xc(m+1)
            y(j)=3*(((x(j)-xc(m))/(xc(m+1)-xc(m)))^2)-...
                2*(((x(j)-xc(m))/(xc(m+1)-xc(m)))^3);
        end
    else
        if xc(l)<x(j) && x(j)<=xc(l+1)
            y(j)=3*(((x(j)-xc(l))/(xc(l+1)-xc(l)))^2)-...
                2*(((x(j)-xc(l))/(xc(l+1)-xc(l)))^3);
        elseif xc(l+1)<=x(j) && x(j)<xc(l+2)
            y(j)=3*(((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)))^2)-...
                2*(((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)))^3);
        end
    end
end
end
function y=fS(x,l,xc) %non toccare
m=length(xc)-1;
y=zeros(size(x));
for j=1:length(x)
    if l==0
        if xc(1)<=x(j) && x(j)<xc(2)
            y(j)=((xc(2)-x(j))^2/(xc(2)-xc(1)))-...
                ((xc(2)-x(j))^3/(xc(2)-xc(1))^2);
        end
    elseif l==m
        if xc(m)<x(j) && x(j)<=xc(m+1)
            y(j)=-((x(j)-xc(m))^2/(xc(m+1)-xc(m)))+...
                ((x(j)-xc(m))^3/(xc(m+1)-xc(m))^2);
        end
    else
        if xc(l)<x(j) && x(j)<=xc(l+1)
            y(j)=-((x(j)-xc(l))^2/(xc(l+1)-xc(l)))+...
                ((x(j)-xc(l))^3/(xc(l+1)-xc(l))^2);
        elseif xc(l+1)<=x(j) && x(j)<xc(l+2)
            y(j)=((xc(l+2)-x(j))^2/(xc(l+2)-xc(l+1)))-...
                ((xc(l+2)-x(j))^3/(xc(l+2)-xc(l+1))^2);
        end
    end
end
end
function y=dH(x,l,xc) %non toccare
m=length(xc)-1;
y=zeros(size(x));
for j=1:length(x)
    if l==0
        if xc(1)<=x(j) && x(j)<xc(2)
            y(j)=(6/(xc(2)-xc(1)))*((xc(2)-x(j))/(xc(2)-xc(1)))*...
                ((xc(2)-x(j))/(xc(2)-xc(1))-1);
        end
    elseif l==m
        if xc(m)<x(j) && x(j)<=xc(m+1)
            y(j)=(6/(xc(m+1)-xc(m)))*((x(j)-xc(m))/(xc(m+1)-xc(m)))*...
                (1-((x(j)-xc(m))/(xc(m+1)-xc(m))));
        end
    else
        if xc(l)<x(j) && x(j)<=xc(l+1)
            y(j)=(6/(xc(l+1)-xc(l)))*((x(j)-xc(l))/(xc(l+1)-xc(l)))*...
                (1-((x(j)-xc(l))/(xc(l+1)-xc(l))));
        elseif xc(l+1)<=x(j) && x(j)<xc(l+2)
            y(j)=(6/(xc(l+2)-xc(l+1)))*((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)))*...
                ((xc(l+2)-x(j))/(xc(l+2)-xc(l+1))-1);
        end
    end
end
end
function y=dS(x,l,xc) %non toccare
m=length(xc)-1;
y=zeros(size(x));
for j=1:length(x)
    if l==0
        if xc(1)<=x(j) && x(j)<xc(2)
            y(j)=3*(((xc(2)-x(j))/(xc(2)-xc(1)))^2)-...
                2*((xc(2)-x(j))/(xc(2)-xc(1)));
        end
    elseif l==m
        if xc(m)<x(j) && x(j)<=xc(m+1)
            y(j)=3*(((x(j)-xc(m))/(xc(m+1)-xc(m)))^2)-...
                2*((x(j)-xc(m))/(xc(m+1)-xc(m)));
        end
    else
        if xc(l)<x(j) && x(j)<=xc(l+1)
            y(j)=3*(((x(j)-xc(l))/(xc(l+1)-xc(l)))^2)-...
                2*((x(j)-xc(l))/(xc(l+1)-xc(l)));
        elseif xc(l+1)<=x(j) && x(j)<xc(l+2)
            y(j)=3*(((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)))^2)-...
                2*((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)));
        end
    end
end
end
function y=fu(x,l,xc) %scrivi quello che hai trovato nei calcoli della base u_l, attento all'indice l
    if l==2*m
        y=fS(x,m,xc);
    elseif mod(l,2)==1
        y=fS(x,(l-1)/2,xc);
    else
        y=fH(x,l/2,xc);
    end
end
function y=du(x,l,xc) %derivata di u_l
    if l==2*m
        y=dS(x,m,xc);
    elseif mod(l,2)==1
        y=dS(x,(l-1)/2,xc);
    else
        y=dH(x,l/2,xc);
    end
end
function [lb,rb]=support(l,xc)
    %se abbiamo u_0 aggiungo if l==0 lb=xc(1); rb=xc(2);
    if l>=2*m 
        lb=xc(m);
        rb=xc(m+1);
    elseif mod(l,2)==1 %sono le funzioni S
        %In generale alcuni supporti sono più piccoli ma questo va bene per
        %tutti, idem nelle funzioni H
        lb=xc((l-1)/2);
        rb=xc(2+(l-1)/2);
    else %sono le funzioni H
        lb=xc(l/2);
        rb=xc(2+l/2);
    end
end
function y=fA(x,l,j,xc)
%Inserisci la funzione integranda per la matrice A
y=(x.^2).*du(x,l,xc).*du(x,j,xc)+20.*fu(x,l,xc).*fu(x,j,xc);
end
function y=fb(x,l,xc)
%Inserisci la funzione integranda per il vettore b
y=(7.*(x.^2)-5).*fu(x,l,xc);
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2023-11-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2043132-energy-method-with-hermite-cubic-function-2-essential-condition#answer_1346867

m and nv are inputs to the function "ex_enerHermite".

Where do you set these values and where do you call the function ?

5 个评论
显示 3更早的评论隐藏 3更早的评论

Aurora 2023-11-5

在 MATLAB Online 中打开

Hello, the problem is change. The approximatly solution is a littlet different.

In the followinf, you can see the listing and the output.

function [x,u]=ex_enerHermite2(m,nv)
%trasla il problema se necessario con condizione essenziale in b
B=11/14;
C=-6/7;
h=1.5/m;
xc=-2:h:-0.5;
A=zeros(2*m+1);%se non ho condizioni essenziali diventa 2*m+2
b=zeros(2*m+1,1);%se non ho condizioni essenziali diventa 2*m+2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Caso condizione essenziale in a
for l=1:2*m+1 %se non ho condizioni essenziali diventa 2*m+2
    [lbl,rbl]=support(l,xc);
    for j=1:2*m+1 %se non ho condizioni essenziali diventa 2*m+2
        if abs(l-j)<4 %altrimenti l'integrale Ã¨ zero causa supporto base
            [lbj,rbj]=support(j,xc);
            lb=max(lbl,lbj);
            rb=min(rbl,rbj);
            A(l,j)=integral(@(x)fA(x,l,j,xc),lb,rb);
            
        end
    end
    b(l)=integral(@(x)fb(x,l,xc),lbl,rbl);
end
A(1,1)=2+A(1,1);
%H_m contiene anche il termine fuori dall'integrale
%da aggiungere.
delta=A\b;
h=1.5/nv;
x=-2:h:-0.5;
u=zeros(1,nv+1);
for l=1:2*m+1 %da cambiare in base a dove fai la sommatoria della soluzione u
    u=u+delta(l)*fu(x,l,xc);
    %se cond. ess. in a allora l=1:2*m+1 e 
    %u=u+delta(l)*fu(x,l,xc)  
    %se non ho cond. ess. allora l=0:2*m+1 e 
    %u=u+delta(l+1)*fu(x,l,xc) 
end
%u Ã¨ la soluzione del problema traslato quindi se necessario aggiungo l(x)
u=u+B*x+C;
end
function y=fH(x,l,xc) %non toccare
m=length(xc)-1;
y=zeros(size(x));
for j=1:length(x)
    if l==0
        if xc(1)<=x(j) && x(j)<xc(2)
            y(j)=3*(((xc(2)-x(j))/(xc(2)-xc(1)))^2)-...
                2*(((xc(2)-x(j))/(xc(2)-xc(1)))^3);
        end
    elseif l==m
        if xc(m)<x(j) && x(j)<=xc(m+1)
            y(j)=3*(((x(j)-xc(m))/(xc(m+1)-xc(m)))^2)-...
                2*(((x(j)-xc(m))/(xc(m+1)-xc(m)))^3);
        end
    else
        if xc(l)<x(j) && x(j)<=xc(l+1)
            y(j)=3*(((x(j)-xc(l))/(xc(l+1)-xc(l)))^2)-...
                2*(((x(j)-xc(l))/(xc(l+1)-xc(l)))^3);
        elseif xc(l+1)<=x(j) && x(j)<xc(l+2)
            y(j)=3*(((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)))^2)-...
                2*(((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)))^3);
        end
    end
end
end
function y=fS(x,l,xc) %non toccare
m=length(xc)-1;
y=zeros(size(x));
for j=1:length(x)
    if l==0
        if xc(1)<=x(j) && x(j)<xc(2)
            y(j)=((xc(2)-x(j))^2/(xc(2)-xc(1)))-...
                ((xc(2)-x(j))^3/(xc(2)-xc(1))^2);
        end
    elseif l==m
        if xc(m)<x(j) && x(j)<=xc(m+1)
            y(j)=-((x(j)-xc(m))^2/(xc(m+1)-xc(m)))+...
                ((x(j)-xc(m))^3/(xc(m+1)-xc(m))^2);
        end
    else
        if xc(l)<x(j) && x(j)<=xc(l+1)
            y(j)=-((x(j)-xc(l))^2/(xc(l+1)-xc(l)))+...
                ((x(j)-xc(l))^3/(xc(l+1)-xc(l))^2);
        elseif xc(l+1)<=x(j) && x(j)<xc(l+2)
            y(j)=((xc(l+2)-x(j))^2/(xc(l+2)-xc(l+1)))-...
                ((xc(l+2)-x(j))^3/(xc(l+2)-xc(l+1))^2);
        end
    end
end
end
function y=dH(x,l,xc) %non toccare
m=length(xc)-1;
y=zeros(size(x));
for j=1:length(x)
    if l==0
        if xc(1)<=x(j) && x(j)<xc(2)
            y(j)=(6/(xc(2)-xc(1)))*((xc(2)-x(j))/(xc(2)-xc(1)))*...
                ((xc(2)-x(j))/(xc(2)-xc(1))-1);
        end
    elseif l==m
        if xc(m)<x(j) && x(j)<=xc(m+1)
            y(j)=(6/(xc(m+1)-xc(m)))*((x(j)-xc(m))/(xc(m+1)-xc(m)))*...
                (1-((x(j)-xc(m))/(xc(m+1)-xc(m))));
        end
    else
        if xc(l)<x(j) && x(j)<=xc(l+1)
            y(j)=(6/(xc(l+1)-xc(l)))*((x(j)-xc(l))/(xc(l+1)-xc(l)))*...
                (1-((x(j)-xc(l))/(xc(l+1)-xc(l))));
        elseif xc(l+1)<=x(j) && x(j)<xc(l+2)
            y(j)=(6/(xc(l+2)-xc(l+1)))*((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)))*...
                ((xc(l+2)-x(j))/(xc(l+2)-xc(l+1))-1);
        end
    end
end
end
function y=dS(x,l,xc) %non toccare
m=length(xc)-1;
y=zeros(size(x));
for j=1:length(x)
    if l==0
        if xc(1)<=x(j) && x(j)<xc(2)
            y(j)=3*(((xc(2)-x(j))/(xc(2)-xc(1)))^2)-...
                2*((xc(2)-x(j))/(xc(2)-xc(1)));
        end
    elseif l==m
        if xc(m)<x(j) && x(j)<=xc(m+1)
            y(j)=3*(((x(j)-xc(m))/(xc(m+1)-xc(m)))^2)-...
                2*((x(j)-xc(m))/(xc(m+1)-xc(m)));
        end
    else
        if xc(l)<x(j) && x(j)<=xc(l+1)
            y(j)=3*(((x(j)-xc(l))/(xc(l+1)-xc(l)))^2)-...
                2*((x(j)-xc(l))/(xc(l+1)-xc(l)));
        elseif xc(l+1)<=x(j) && x(j)<xc(l+2)
            y(j)=3*(((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)))^2)-...
                2*((xc(l+2)-x(j))/(xc(l+2)-xc(l+1)));
        end
    end
end
end
function y=fu(x,l,xc) %scrivi quello che hai trovato nei calcoli della base u_l, attento all'indice l
m = length(xc)-1;  
     if l==2*m+1
        y=fS(x,m,xc);
     elseif mod(l,2)==1
        y=fS(x,(l-1)/2,xc);
    else
        y=fH(x,(l-2)/2,xc);
    end
end
function y=du(x,l,xc) %derivata di u_l
m = length(xc)-1;  
     if l==2*m+1
         y=dS(x,m,xc);
    elseif mod(l,2)==1
        y=dS(x,(l-1)/2,xc);
    else
        y=dH(x,(l-2)/2,xc);
    end
end
function [lb,rb]=support(l,xc)
    m=length(xc)-1;
    %se abbiamo u_0 aggiungo if l==0 lb=xc(1); rb=xc(2);
    if l==1 || l==2
        lb=xc(1);
        rb=xc(2);
    elseif l==2*m+1 %cioè il fondo, se non ho u_2m+1 posso scrivere if l==2*m
        lb=xc(m);
        rb=xc(m+1);
    elseif mod(l,2)==1 %sono le funzioni S
        %In generale alcuni supporti sono più piccoli ma questo va bene per
        %tutti, idem nelle funzioni H
        lb=xc((l-1)/2);
        rb=xc(2+(l-1)/2);
    else %sono le funzioni H
        lb=xc((l-2)/2);
        rb=xc(2+(l-2)/2);
    end
end
function y=fA(x,l,j,xc)
%Inserisci la funzione integranda per la matrice A
y=(x.^2).*du(x,l,xc).*du(x,j,xc)+6.*fu(x,l,xc).*fu(x,j,xc);
end
function y=fb(x,l,xc)
B=11/14;
C=-6/7;
%Inserisci la funzione integranda per il vettore b
y=(8.*x+1./x-4*B.*x-6*C).*fu(x,l,xc);
end

Aurora 2023-11-5

You can choose m, nv like you prefer. In particoular, i choose 15,6

Aurora 2023-11-5

在 MATLAB Online 中打开

>> [x,u]=ex_enerHermite2(30,6)
x =
   -2.0000   -1.7500   -1.5000   -1.2500   -1.0000   -0.7500   -0.5000
u =
   -2.9103   -2.6970   -2.4510   -2.1745   -1.8722   -1.5543   -1.2500
>> y=((x.^2)/3)+(2.*x)+1./(6.*x)
y =
   -2.7500   -2.5744   -2.3611   -2.1125   -1.8333   -1.5347   -1.2500
>> plot(x,u,'r*',x,y)