Element of order 2

Construction and an example of an element of order 2

Define the interpolated function

Clear[f, a, b, c, x]\(f[x_]\  = \ a*x^2\  + \ b*x\  + c; 

    fInterpol[x_] = 
      f[x] /. 
		
        Solve[{f[\(-h\)] == yl, f[0] == yc, f[h] == yr}, {a, b, c}]; 

    Needs["NumericalMath`GaussianQuadrature`"]\)\(\({{pl, pc, pr}, {wl, wc, wr}} = 
	
      Transpose[GaussianQuadratureWeights[3, \ \(-h\), \ h]]; \)\)

Interpolation of the function

Evaluate the function at the Gauss quadrature points and represent the formula using the notation of a matrix.

\(m = Simplify[Map[fInterpol, {pl, pc, pr}]]; 

    M0 = Transpose[Map[Flatten[Coefficient[m, #]]&, {yl, yc, yr}]]; \)

Test the code and compare with the lecture notes.

M0 . {yl, yc, yr}

Transpose[M0] . {5, 8, 5}/9/2

Interpolation of the derivative of function

Evaluate the function at the Gauss quadrature points and represent the formula using the notation of a matrix.

\(m = Simplify[Map[fInterpol', {pl, pc, pr}]]; 

    M1 = Transpose[Map[Flatten[Coefficient[m, #]]&, {yl, yc, yr}]]; \)

Test the code and compare with the lecture notes.

M1 . {yl, yc, yr}

Construction of the elements stiffnes matrix

\(ElementStiffnessMatrix[a_, b_, {xl_, xr_}]\  := 
	
      Module[{pl, pc, pr, w, h1, M1h, M0h, avalues, bvalues}, 

        		{{pl, pc, pr}, w} = 
	
          Transpose[GaussianQuadratureWeights[3, \ xl, \ xr]]; 
		
        h1 = \((xr - xl)\)/2; 
		M0h\  = \ M0 /. \ h -> h1; 
		
        M1h\  = \ M1 /. \ h -> h1; 
		avalues = Map[a, {pl, pc, pr}]; 

        		bvalues = Map[b, {pl, pc, pr}]; 
		
        Transpose[M1h] . DiagonalMatrix[avalues*{5, 8, 5}/9*h1] . M1h + 
		
          Transpose[M0h] . DiagonalMatrix[bvalues*{5, 8, 5}/9*h1] . M0h
	]
        \)\(ElementRHS[f_, {xl_, xr_}]\  := 
	
      Module[{pl, pc, pr, w, h1, M0h, fvalues}, 

        		{{pl, pc, pr}, w} = 
	
          Transpose[GaussianQuadratureWeights[3, \ xl, \ xr]]; 
		
        h1 = \((xr - xl)\)/2; 
		M0h\  = \ M0 /. \ h -> h1; 
		
        fvalues = Map[f, {pl, pc, pr}]; 
		
        Transpose[M0h] . fvalues*{5, 8, 5}/9*h1
	]\)

A simple test case

\(a[x_] = 1; 
b[x_] = 0.1; 
f[x_] = 1; 

    ElementStiffnessMatrix[a, b, {\(-3\), 2}] // MatrixForm\)ElementRHS[f, {-3, 2}] // MatrixForm

Exact solution of the test problem

Clear[A, x, F, Em, eps]A[x_]  :=  10 - 0.09 x\(F = 2*10^4; 
Em = 3*10^6; 
eps[x_] = F/\((Em*A[x])\)\)uexact[x_] =  Integrate[eps[s], {s, 0, x}]

uexact[100]

Piecewise quadratic interpolation code

\(interpol2[x_, data_] := \ 
      Module[{xp, yp, pos, f}, 
		{xp, yp} = Transpose[data]; 
	
        pos = Floor[\((Length[Select[xp, \((# <= x)\)&]] + 1)\)/2]; 
	
        f = Interpolation[data[\([{2*pos - 1, 2*pos, 2*pos + 1}]\)], 
            InterpolationOrder\  -> \ 2]; 
	f[x]]\)\(Dinterpol2[x_, data_] := \ 
      Module[{xp, yp, pos, f}, 
		{xp, yp} = Transpose[data]; 
	
        pos = Floor[\((Length[Select[xp, \((# <= x)\)&]] + 1)\)/2]; 
	
        f = Interpolation[data[\([{2*pos - 1, 2*pos, 2*pos + 1}]\)], 
            InterpolationOrder\  -> \ 2]; 
	\(f'\)[x]]\)

A simple test of the above routines

\(Plot[interpol2[x, {{0, 0}, {1, 1}, {2, 0}, {3, 1}, {4, 0}}], 

      	{x, 0, 4}]; 

    Plot[Dinterpol2[x, {{0, 0}, {1, 1}, {2, 0}, {3, 1}, {4, 0}}], 

      	{x, 0, 4}]; \)

The test problem solved with elements of order 2

Solution with 2 elements

Clear[A, b, x, x0]\(x0 = 0; x2 = 50; x4 = 100; 
x = {0, 25, 50, 75, 100}; 

    A[xt_] = \((10 - 0.09*xt)\)*3*10^6; 
b[xt_] = 0; 

    M = Table[0, {5}, {5}]; 


    mElement = N[ElementStiffnessMatrix[A, b, {x0, x2}]]; 

    For[i = 1, i <= 3, \(i++\), 
	\ 
      For[j = 1, j <= 3, \(j++\), 
		
        M[\([i, j]\)] += mElement[\([i, j]\)]]]; 


    mElement = N[ElementStiffnessMatrix[A, b, {x2, x4}]]; 

    For[i = 1, i <= 3, \(i++\), 
	\ 
      For[j = 1, j <= 3, \(j++\), 
		
        M[\([2 + i, 2 + j]\)] += mElement[\([i, j]\)]]]; 


    RHS = Table[0, {5}]; RHS[\([5]\)] = 2*10^4; 

    M[\([1]\)] = {1, 0, 0, 0, 0}; 
u = LinearSolve[M, RHS]; 


    data = Transpose[{x, u}]; \)

\(\(Plot[{uexact[xt], interpol2[xt, data]}, {xt, 0, 100}, 
      PlotRange\  -> All, 	
		AxesLabel -> {"x", "Strain"}, 

      PlotStyle -> {Dashing[{0.03, 0.02}], Dashing[{}]}]; \)\)

\(Duexact[xt_] = D[uexact[xt], xt]; 

    Plot[{Duexact[xt], Dinterpol2[xt, data]}, {xt, 0, 100}, 
      PlotRange\  -> {0, 0.007}, 	
		
      AxesLabel -> {"x", "Strain"}, 

      PlotStyle -> {Dashing[{0.03, 0.02}], Dashing[{}]}]; \)

Solution with 5 elements

\(x = Table[k, {k, 0, 100, 10}]; 

    A = Function[x, \((10 - 0.09*x)\)*3*10^6]; 
M = Table[0, {11}, {11}]; 

    For[e = 1, e <= 5, \(e++\), 

      mElement = N[ElementStiffnessMatrix[A, b, {\((2*e - 2)\)*10, 2*e*10}]]; 
      
For[i = 1, i <= 3, \(i++\), 
	\ 
        For[j = 1, j <= 3, \(j++\), 
		
          M[\([\((e - 1)\)*2 + i, \((e - 1)\)*2 + j]\)] += 
            mElement[\([i, j]\)]]]; ]
\)\(RHS = Table[0, {11}]; 
RHS[\([11]\)] = 2*10^4; 

    M[\([1]\)] = Table[0, {11}]; 
M[\([1, 1]\)] = 1; 

    u = LinearSolve[M, RHS]; 

data = Transpose[{x, u}]; \)

\(\(Plot[{uexact[xt], interpol2[xt, data]}, {xt, 0, 100}, 
      PlotRange\  -> All, 	
		AxesLabel -> {"x", "Displacement"}, 
      
PlotStyle -> {Dashing[{0.03, 0.02}], Dashing[{}]}]; \)\)

\(\(Plot[
      \((uexact[xt] - interpol2[xt, data])\)/uexact[100], {xt, 0, 100}, 

      PlotRange\  -> All]; \)\)

\(Duexact[xt_] = D[uexact[xt], xt]; 

    Plot[{Duexact[xt], Dinterpol2[xt, data]}, {xt, 0, 100}, 
      PlotRange\  -> {0, 0.007}, 	
		
      AxesLabel -> {"x", "Strain"}, 

      PlotStyle -> {Dashing[{0.03, 0.02}], Dashing[{}]}]; \)

Stand alone code

The code to be included

\(\(\( (*\(SetDirectory["~sha/tex/wahlfach/Variationsrechnung"];\)*) \)\(
\)\(<< BVP.m\)\)\)

Test with BC

?BVP

Test1

<< BVP.m

Clear[a, b, c, f, x, n]
a = Function[x, -1];
b = Function[x, 2];
f = Function[x, -x^3];
n = 3;
x = Table[t, {t, 0, 1, 1/n}];
data = BVP[a, b, f, x, {"D", "D"}, {0, 0.1}];
\(\(g1 = 
        ListPlot[
          data, 
		PlotStyle\  -> \ 
            PointSize[0.02], 
		DisplayFunction\  -> 
            Identity]\ ;\)\)
\(\(yn[t_] = 
        y[t] /. DSolve[
			{D[a[t]*\(y'\)[t], t] - b[t]*y[t] == 
                f[t], 
				y[0] == 0, y[1] == 0.1}, y[t], 
            t];\)\)
g2 = 
        Plot[yn[t], {t, 0, 1}, DisplayFunction  ->  Identity];
\(\(Show[{g1, g2}, 
	PlotRange\  -> \ All, 
	
        DisplayFunction\  -> $DisplayFunction];\)\)

<< Interpol2.m
Clear[x]
yfem[x_] := Interpol2[x, data];
gy = Plot[{yfem[x], yn[x]}, {x, 0, 1}];
dyfem[x_] := DInterpol2[x, data];
dyn[x_] = D[yn[x], x];
\(\(gdy = 
        Plot[{dyfem[x], dyn[x]}, {x, 0, 1}, 

          PlotRange\  -> All];\)\)

Test2

Clear[a, b, c, f, x, t, n]
a = Function[x, x];
b = Function[x, 0];
\(\(f = 
        Function[x, Which[x < 1, 0, x < 2, \(-x\), True, 0]];\)
 (*\(f = 
          Function[x, 
            x*\((Sign[x - 1] + 1)\)*\((Sign[x - 2] - 1)\)/4];\)*) \)
n = 15;
x = Table[t, {t, 0, 3, 3/n}];
data = BVP[a, b, f, x, {"N", "D"}, {0, 0}];
\(\(g1 = 
        ListPlot[
          data, 
		PlotStyle\  -> \ 
            PointSize[0.02], 
		PlotJoined\  -> \ True];\)\)

Test3

Clear[x]
sol[x_]  =  (1 - x^2)^2;
ff[x_] = 
        Simplify[D[ (x + 1)*sol'[x], x] +  sol[x]];

Clear[a, b, c, f, x, t, n]
a = Function[x, 1 + x];
b = Function[x, -1];
f = Function[x, ff[x]];
n = 2;
x = Table[t, {t, 0, 1, 1/n}];
data = BVP[a, b, f, x, {"N", "N"}, {0, 0}];
\(\(g1 = 
        ListPlot[
          data, 
		PlotStyle\  -> \ 
            PointSize[0.02], 
		DisplayFunction\  -> 
            Identity]\ ;\)\)
\(	\(yn[t_] = 
        y[t] /. NDSolve[
			{D[a[t]*\(y'\)[t], t] - b[t]*y[t] == 
                f[t], 
				\(y'\)[0] == 0, \(y'\)[1] == 0}, 
            y[t], {t, 0, 1}];\)\)
g2 = 
        Plot[{sol[t], yn[t]}, {t, 0, 1}, 
          DisplayFunction  ->  Identity];
\(\(Show[{g1, g2}, 
	PlotRange\  -> \ All, 
	
        DisplayFunction\  -> $DisplayFunction];\)\)

Test4

Clear[a, b, c, f, x, n]
a = Function[x, 1 + x];
b = Function[x, -2];
f = Function[x, Tanh[x^2]];
n = 50;
x = Table[t, {t, 0, 4, 4/n}];
data = BVP[a, b, f, x, {"D", "N"}, {17, 0}];
\(\(g1 = 
        ListPlot[
          data, 
		PlotStyle\  -> \ 
            PointSize[0.01], 
		DisplayFunction\  -> 
            Identity]\ ;\)\)
\(\(yn[t_] = 
        y[t] /. NDSolve[
			{D[a[t]*\(y'\)[t], t] - b[t]*y[t] == 
                f[t], 
				y[0] == 17, \(y'\)[4] == 0}, 
            y[t], {t, 0, 4}];\)\)
g2 = 
        Plot[yn[t], {t, 0, 4}, DisplayFunction  ->  Identity];
\(\(Show[{g1, g2}, 
	PlotRange\  -> \ All, 
	
        DisplayFunction\  -> $DisplayFunction];\)\)

Test5

<< BVP.m

Clear[a, b, c, f, x, n]
a = Function[x, -1];
b = Function[x, 2];
f = Function[x, -x^3];
n = 3;
x = Table[t, {t, 0, 1, 1/n}];
data = 
        BVP[a, b, f, x, {"N", "N"}, {+1, +1}];
\(\(g1 = 
        ListPlot[
          data, 
		PlotStyle\  -> \ 
            PointSize[0.02], 
		DisplayFunction\  -> 
            Identity]\ ;\)\)
\(\(yn[t_] = 
        y[t] /. DSolve[
			{D[a[t]*\(y'\)[t], t] - b[t]*y[t] == 
                f[t], 
				\(y'\)[0] == \(-1\), \(y'\)[1] == 1}, y[t], 
            t];\)\)
g2 = 
        Plot[yn[t], {t, 0, 1}, DisplayFunction  ->  Identity];
\(\(Show[{g1, g2}, 
	PlotRange\  -> \ All, 
	
        DisplayFunction\  -> $DisplayFunction];\)\)

<< Interpol2.m
Clear[x]
yfem[x_] := Interpol2[x, data];
gy = Plot[{yfem[x], yn[x]}, {x, 0, 1}];
dyfem[x_] := DInterpol2[x, data];
dyn[x_] = D[yn[x], x];
\(\(gdy = 
        Plot[{dyfem[x], dyn[x]}, {x, 0, 1}, 

          PlotRange\  -> All];\)\)

Test6 (Hausaufgabe)

<< BVP.m

Clear[a, b, c, f, x, n, yn, y, t]
a = Function[x, x];
b = Function[x, -7];
f = Function[x, -Sin[x]];
n = 10;
x = Table[t, {t, 0, 3, 1/n}];
data = BVP[a, b, f, x, {"D", "N"}, {3, 0}];
\(\(g1 = 
        ListPlot[
          data, 
		PlotStyle\  -> \ 
            PointSize[0.01], 
		DisplayFunction\  -> 
            Identity]\ ;\)\)
\(\(yn[t_] = 
        y[t] /. NDSolve[
			{D[a[t]*\(y'\)[t], t] - b[t]*y[t] == 
                f[t], 
				y[data[\([2, 1]\)]] == 
                data[\([2, 2]\)], \(y'\)[3] == 0}, 
            y[t], {t, 0.01, 3}];\)\)
g2 = 
        Plot[yn[t], {t, data[[2, 1]], 1}, 
          DisplayFunction  ->  Identity];
\(\(Show[{g1, g2}, 
	PlotRange\  -> \ All, 
	
        DisplayFunction\  -> $DisplayFunction];\)\)

<< Interpol2.m
Clear[x]
yfem[x_] := Interpol2[x, data];
gy = Plot[{yfem[x], yn[x]}, {x, data[[2, 1]], 1}];
dyfem[x_] := DInterpol2[x, data];
dyn[x_] = D[yn[x], x];
\(\(gdy = 
        Plot[{dyfem[x], dyn[x]}, {x, data[\([2, 1]\)], 1}, 

          PlotRange\  -> All];\)\)

The test problem solved with elements of order 2

<< BVP.m
<< Interpol2.m

Clear[A, b, f, x, data]
A = Function[x, -(10 - 0.09*x)*3*10^6];
b = Function[x, 0];
f = Function[x, 0];
n = 5;
uexact[
        x_] =  -Integrate[
          2*10^4/A[s], {s, 0, x}]
x = Table[s, {s, 0, 100, 100/n}];
data = 
        BVP[A, b, f, x, {"D", "N"}, {0, 2*10^4}];

\(\(g1 = 
        Plot[{uexact[xt], Interpol2[xt, data]}, {xt, 0, 100}, 
          PlotRange\  -> 
            All, 	
		AxesLabel -> {"x", "Displacement"}, 

          PlotStyle -> {Dashing[{0.03, 0.02}], Dashing[{}]}];\)\)

Duexact[xt_] = D[uexact[xt], xt];
\(\(g2 = 
        Plot[{Duexact[xt], DInterpol2[xt, data]}, {xt, 0, 100}, 
          PlotRange\  -> {0, 
              0.007}, 	
		AxesLabel -> {"x", "Strain"}, 

          PlotStyle -> {Dashing[{0.03, 0.02}], Dashing[{}]}];\)\)

\(\(g3 = 
        Plot[\((uexact[xt] - Interpol2[xt, data])\)/uexact[100], {xt, 0, 
            100}, 
PlotRange\  -> All];\)\)

Convergence Test 1

<< BVP.m;
<< Interpol2.m;

Clear[a, b, c, f, x, t, n];
a = Function[x, 2];
b = Function[x, -1];
f = Function[x, Sin[Pi/2 - x]];
n = 22;
x = Table[t, {t, 0, 2, 2/n}];
data = BVP[a, b, f, x, {"N", "D"}, {0, 0}];
\(\(g1 = 
        ListPlot[data, 
		PlotStyle\  -> \ PointSize[0.01]];\)\)
Clear[x]
\(\(uexact[x_] = 
        y[x] /. \(DSolve[{D[a[x]*D[y[x], x], x] - b[x]*y[x] == 
                  f[x], \(y'\)[0] == 0, y[2] == 0}, 
y[x], 
              x]\)[\([1]\)];\)
	\)
Duexact[x_] = D[uexact[x], x];

un[x_]  :=  Interpol2[x, data]
Plot[{un[x], uexact[x]}, {x, 0, 2}, PlotRange -> All];

testx =  1.01;
\(\(res = 
        Table[{n = 3*k + 1, 
            x = Table[t, {t, 0, 2, 2/n}]; 
\ \ \ \ \ \ \ \ \ \ \ \ data = 
              BVP[a, b, f, 
                x, {"N", "D"}, {0, 0}]; 
			\ \ \ \ \ \ Abs[
              Interpol2[testx, data] - uexact[testx]], 
				Abs[
              DInterpol2[testx, data] - Duexact[testx]]}, 
			{k, 1, 30, 
            3}];\)\)
\(\(resu = 
        Map[Drop[#, \(-1\)] &, res];\)
 (*\(ListPlot[resu, 
          PlotJoined\  -> True, PlotRange\  -> All];\)*) \)
ListPlot[Log[resu], PlotJoined  -> True, 
        PlotRange -> All];
logresu = Log[Abs[resu]];
n = Length[resu];
(logresu[[1, 2]] - 
          logresu[[n, 2]])/(logresu[[1, 1]] - 
          logresu[[n, 1]])

testx =  1.01;
\(\(res = 
        Table[{n = 3*k + 1, 
            x = Table[t, {t, 0, 2, 2/n}]; 
\ \ \ \ \ \ \ \ \ \ \ \ data = 
              BVP[a, b, f, 
                x, {"N", "D"}, {0, 0}]; 
			\ \ \ \ \ \ Abs[
              Interpol2[testx, data] - uexact[testx]], 
				Abs[
              DInterpol2[testx, data] - Duexact[testx]]}, 
			{k, 1, 30, 
            3}];\)\)
\(\(resu = 
        Map[Drop[#, \(-1\)] &, res];\)
 (*\(ListPlot[resu, 
          PlotJoined\  -> True, PlotRange\  -> All];\)*) \)
g2 = 
        ListPlot[Log[resu], PlotJoined  -> True, 
          PlotRange -> All];
logresu = Log[Abs[resu]];
n = Length[resu];
(logresu[[1, 2]] - 
          logresu[[n, 2]])/(logresu[[1, 1]] - 
          logresu[[n, 1]])

\(\(resu = 
        Map[Drop[#, {2}] &, res];\)
 (*\(ListPlot[resu, 
          PlotJoined\  -> True, PlotRange\  -> All];\)*) \)
ListPlot[Log[resu], PlotJoined  -> True, 
        PlotRange -> All];

logresu = Log[Abs[resu]];
n = Length[resu];
(logresu[[1, 2]] - 
          logresu[[n, 2]])/(logresu[[1, 1]] - 
          logresu[[n, 1]])

Convergence Test 2

<< BVP.m;
<< Interpol2.m;

Clear[a, b, c, f, x, t, n];
a = Function[x, 2];
b = Function[x, -1];
f = Function[x, -x^2 + x^3];
n = 22;
x = Table[t, {t, 0, 2, 2/n}];
data = BVP[a, b, f, x, {"D", "D"}, {0, 0}];
\(\(g1 = 
        ListPlot[data, 
		PlotStyle\  -> \ PointSize[0.01]];\)\)
Clear[x]
\(uexact[x_] = 
      y[x] /. \(DSolve[{D[a[x]*D[y[x], x], x] - b[x]*y[x] == f[x], y[0] == 0, 
              y[2] == 0}, 
y[x], x]\)[\([1]\)]\)

un[x_]  :=  Interpol2[x, data]
Plot[{un[x], uexact[x]}, {x, 0, 2}, PlotRange -> All];

testx = 1.0;
\(res = 
      Table[{n = 3*k + 1, 
          x = Table[t, {t, 0, 2, 2/n}]; 
\ \ \ \ \ \ \ \ \ \ \ \ data = 
            BVP[a, b, f, 
              x, {"D", "D"}, {0, 0}]; 
			\ \ \ \ \ \ Abs[
            Interpol2[testx, data] - uexact[testx]]}, 
			{k, 1, 30, 
          3}]\)

ListPlot[res, PlotJoined  -> True, PlotRange  -> All];

ListPlot[Log[res], PlotJoined  -> True, 
        PlotRange -> All];

logres = Log[Abs[res]];
n = Length[res];
(logres[[1, 2]] - logres[[n, 2]])/(logres[[1, 1]] - 
          logres[[n, 1]])

Construction of an element of order 4

Define the interpolated function

Clear[f, x]
f[x_]  =  a0  +  a1*x + a2*x^2 + a3*x^3 + a4*x^4;
\(\(fInterpol[x_] = 
        f[x] /. 
		Solve[{f[\(-h\)] == yl, \(f'\)[\(-h\)] == yl1, 
              f[0] == yc, f[h] == yr, \(f'\)[h] == yr1}, {a0, a1, a2, a3, 
              a4}];\)\)

Interpolation of the function

Evaluate the function at the Gauss quadrature points and represent the formula using the notation of a matrix.

Needs["NumericalMath`GaussianQuadrature`"]
\(\({{pl, pc, pr}, {wl, wc, wr}} = 
	
        Transpose[GaussianQuadratureWeights[3, \ \(-h\), \ h]];\)\)
m = Simplify[Map[fInterpol, {pl, pc, pr}]];
M0 = 
        Transpose[
          Map[Flatten[Coefficient[m, #]] &, {yl, yl1, yc, yr, 
              yr1}]];

Test the code and compare with the lecture notes.

MatrixForm[M0]

Interpolation of the derivative of function

Evaluate the function at the Gauss quadrature points and represent the formula using the notation of a matrix.

m = Simplify[Map[fInterpol', {pl, pc, pr}]];
M1 = 
        Transpose[
          Map[Flatten[Coefficient[m, #]] &, {yl, yl1, yc, yr, 
              yr1}]];

Test the code and compare with the lecture notes.

MatrixForm[M1]

Interpolation of the second derivative of function

Evaluate the function at the Gauss quadrature points and represent the formula using the notation of a matrix.

m = Simplify[Map[fInterpol'', {pl, pc, pr}]];
M2 = 
        Transpose[
          Map[Flatten[Coefficient[m, #]] &, {yl, yl1, yc, yr, 
              yr1}]];

Test the code and compare with the lecture notes.

MatrixForm[M2]

Evaluation of the integrals

K2 := 
        Transpose[M0] . DiagonalMatrix[Map[a2, {pl, pc, pr}]*{wl, wc, wr}] . 
          M0;
K1 := 
        Transpose[M1] . DiagonalMatrix[Map[a1, {pl, pc, pr}]*{wl, wc, wr}] . 
          M1;
K0 := 
      Transpose[M0] . DiagonalMatrix[Map[a0, {pl, pc, pr}]*{wl, wc, wr}] . 
        M0

a2[x_] = 1; a1[x_] = 0; a0[x_] = 0;
h = 1;
(K = K2 + K1 + K0) // MatrixForm;

LinearSolve[K, Transpose[M0] . {0, 1, 0}]

NullSpace[K]

Plot[res[x], {x, -1, 1}]