[资料共享]1D有限元程序Ritz方法，需要解决一个问题，希望探讨

clc; clear; m3bCZ9iE  
N=5; .8(OT./  
xstart=0; xend=4;   %x range (o5j'2:.  
xx=linspace(xstart, xend, N); t.p~\6
Yi  
syms x y F1 F2 ntFT>g{B  
%construct y iOAbaPN  
y=ones(1,N); /D!;u]  
y=sym(y); nKa$1RMO  
for i=1:N ZJPmR/OV_  
    y(i)=['y', num2str(i)] i|`dWOVb  
end J(DN!  
y(1)=0; y(N)=20;     %bound 6;'dUGvH  
f=x^2 + 1/2*x + 1/4; a2IgC25  
for i=1:N-1 f| _u7"OX  
    Fd(i) = 0.5*int( ( ( y(i+1)-y(i) ) / ( xx(i+1)-xx(i) ) )^2, x, xx(i), xx(i+1) )... >P6BW  
    +int(  f * (  y(i) * ( xx(i+1)-x ) / ( xx(i+1)-xx(i) )  + y(i+1) * ( x-xx(i) ) / ( xx(i+1)-xx(i) )   ) , x, xx(i), xx(i+1) );  &i+Ce  
end zk-.u}RBFG  
F=sum(Fd) }}v9 `F  
% how to solve use an equation matrix??? 0:w"M<80  
% F1=[diff(F, y(2)), diff(F, y(3)), diff(F, y(4))] >NYW{(j  
% var1=[y(2), y(3), y(4)] F8-?dpf'  
% var1=solve( F1, var1 ) a=cvCf  
[y(2), y(3), y(4)]=solve( diff(F, y(2)), diff(F, y(3)), diff(F, y(4)), y(2), y(3) , y(4)); ljTBvU  
y H0m|1 7  
for i=1:N `F<jLU^3  
    yy(i)=double(y(i)); Od f[*  
end ;OqB5qd  
plot( xx, yy, 'b*'); hold on ktIi$v  
plot( xx, yy,'b--'); A,  3bC  
%Analytic solution %\]*OZ7  
ax=xstart:0.001:xend; HK`I\,K  
ay=1/12*ax.^4 + 1/12*ax.^3 + 1/8*ax.^2 - 13/6*ax; h8Yx#4  
plot(ax, ay, 'r')
8>hwK)av  
(hiyNMC
  

未注册仅能浏览部分内容，查看全部内容及附件请先 登录 或 注册

感觉楼主的思路有些不对 ]@phF _  
一般来讲 在有限元方法代码中 不需要定义符号变量sym / G7vwC  
需要做的是 你需要知道所有element nodal coordinates 拿1D举例 假设你的computation domain x=(0,1) 用4个element  则每个点的坐标为 0, 0.25, 0.5, 0.75, 1  |'B7v i)  
之后你需要手算出矩阵中每个entry的值(in terms of x) 之后 (jin 的书里面已经给出详细的表达式) 然后将坐标代入即可 'h,VR=e<  
这样子你就会得到刚度矩阵